CC BY
34
данном интервале формируемые прямая /2(а2£ — 0) и обратная ф[ (а1£ + 0) волны представляют волны сжатия.
С увеличением жесткости упругого элемента формируемая в стержне 2 волна деформации /2(а2£ — 0) приближается по форме к падающей волне, уменьшается ее длительность, увеличивается максимальное значение.
На интервале 0 < £ < Т с учетом, что / 1(а1I — 0) = I/Т,
f - о)=r+r ■ ai (/; -о) - r+i
1 - exp (-k1
kT
t
0 < < 1. T
Если жесткость упругого элемента к стремится к бесконечности, то слагаемое
r + 1
1 - exp (-kt
~ 2 г а1 ~
^ стремится к нулю, а функция /2(а2£—0) стремится к значению + ^ — /1 (а1£—0).
Если стержни из одного материала (а1 = а2) и имеют равные волновые сопротивления (г=1), то при к ^ ж функция /2 (а2£ — 0) ^ /1(а1 £ — 0).
Это означает, что при к ^ ж формируемая в стержне 2 волна деформации/(а2£ — 0) стремится воспроизвести падающую волну /1(а1 £ — 0).
Библиографический список
1. Алимов О. Д., Манжосов В. К., Еремьянц В. Э. Распространение волн деформаций в ударных системах. М., 1985. 354 с.
2. Дворников Л. Т., Жуков И. А. Продольный удар по-лукатеноидальным бойком. Новокузнецк, 2006. 80 с.
3. Манжосов В. К. Продольный удар. Ульяновск, 2006. 358 с.
УДК 539.3+622.831
4. Еремьянц В. Э, Невенчанный Ю. В., Писарен-ко Н. Г. Ударное нагружение оснащенных стержней. Фрунзе, 1987. 165 с.
5. Саруев Л. А., Слистин А. П., Авдеева А. И. Передача энергии по ставу штанг при продольном импульсном воздействии. Томск, 1995. 6 с. Деп. в ВИНИТИ 29.11.95, № 3164-В95.
графовый подход при построении
конечно-элементной модели упругих тел
при осесимметричном деформировании
А. А. Тырымов
Волгоградский государственный технический университет E-mail: tyrymov2010@yandex.ru
Предлагается численный метод анализа упругой среды на основе дискретной модели в виде ориентированного графа. В процессе анализа на основе графового подхода тело рассекаем на элементы и для каждого из них строим элементарную ячейку (подграф), являющуюся его моделью. Используя матрицы, представляющие структурные элементы графа, а также уравнения, описывающие разрезанное тело, можно получить уравнения связного тела. Приведены числовые примеры.
Ключевые слова: математическое моделирование, теория упругости, ориентированный граф, деформация, напряжения.
Graph Approach for Finite-Element Based Model of an Elastic Body Under Conditions of Axisymmetric Deformation
A. A. Tyrymov
A numerical method for analysis of the stress -- strain state of elastic media based on a discrete model in form of directed graph is suggested. To analyze a deformable body using the graph approach, we partitione a solid body on elements and replace each element by its model in the form of an elementary cell. The matrices, presenting several structure elements of the graph, and the equations, describing the elementary cells, contribute to deriving the constitutive equations of the intact body. Numerical examples are presented.
Key words: mathematical simulation, elasticity, directed graph, strain, stress.
r
x
Из инженерной практики известно, что элементы многочисленных конструкций имеют форму тел вращения. Во многих случаях геометрия и система нагружения таковы, что исследование напряженного состояния трехмерных тел можно заменить решением соответствующей осесимметричной задачи. В настоящее время одним из самых эффективных и распространенных методов для решения задач механики сплошной среды является метод конечных элементов (МКЭ). Описанию этого метода,
который может быть построен на различных математических идеях, посвящено большое количество монографий. Для вывода основных соотношений МКЭ часто используются основанные на интегральных тождествах методы взвешенных невязок [1], с помощью которых выбором различных весовых функций с единых позиций можно рассмотреть многие, внешне несвязанные, методы (Галеркина, коллокаций, наименьших квадратов, моментов). МКЭ можно сформулировать на основе вариационных принципов минимума полной потенциальной или дополнительной энергии [2]. Одним из самых универсальных в практическом применении для различных сред является принцип виртуальной работы, выраженный в форме виртуальных перемещений. Используются и другие подходы к получению основных зависимостей МКЭ (прямой метод, метод энергетического баланса).
Определение напряженно-деформированного состояния в рамках теории упругости сводится к постановке задачи, включающей в себя помимо граничных условий: дифференциальные уравнения равновесия, которые в случае отсутствия объемных сил имеют вид [Л]т{а} = {0}; геометрические уравнения (формулы Коши) {е} = [Л]{и}; физические соотношения (обобщенный закон Гука) {а} = [С]{е}. Здесь приняты следующие обозначения: [Л] — матрица дифференциальных операторов, связывающая деформации {е} с перемещениями {и}, [С] — матрица модулей упругости материала, {а} — вектор напряжений. Символ Т здесь и далее указывает на операцию транспонирования.
Основополагающее значение при построении МКЭ заключается в специальном выборе системы базисных кусочно-определенных для полученных при дискретизации подобластей Уе функций Ще, каждая из которых имеет конечный носитель. Эти функции должны однозначно восстанавливать некоторую определяющую в данной конкретной задаче функцию в любой точке внутри элемента по дискретным значениям этой функции в узловых точках. Если используется МКЭ в варианте метода перемещений, то за основные неизвестные принимаются компоненты перемещений в узловых точках, расположение которых зависит от формы элемента и вида применяемых для аппроксимации полиномов. Используя определенный набор базисных функций N1, перемещения в произвольной точке элемента представляют в виде {ие} = [Щ]{^е}, где [Щ] — матрица базисных функций конечного элемента, компоненты Щ которой называют функциями формы, {¿е} — вектор узловых перемещений элемента.
С помощью такой аппроксимации и геометрических, а также физических соотношений теории упругости деформации и напряжения в пределах элемента, выражаются через вектор узловых перемещений: {ее} = [Я][Щ{^} = [В]{^}, {ае} = [С][В]{^е}. Здесь [В] = [Л][Щ] — матрица формы деформаций, связывающая деформации с узловыми перемещениями. Варьированием потенциальной энергии конечного элемента или воспользовавшись принципом виртуальной работы [2], устанавливается связь между вектором узловых перемещений и вектором узловых сил, к которому свелись внутренние напряжения {/е} = [Ке]{ие}, где [Ке] = / [В]т[С][В] ¿у — симметричная матрица жесткости конечного элемента.
Следует отметить, что в процессе вывода физические соотношения и геометрические уравнения удовлетворены точно, а уравнения равновесия выполняются в интегральном смысле, т. е. для конечного элемента в целом. При исследовании осесимметричного или плоского деформированного состояния конструкция часто разбивается на четырехугольные конечные элементы. При этом для четырехуз-лового прямоугольного конечного элемента используются билинейные функции формы. Применение «квадратичных» функций формы может быть реализовано только на восьмиузловом прямоугольном элементе [1].
Предлагаемый в работе подход основан на сочетании дискретных и энергетических представлений при моделировании процесса деформирования упругих сред. Дискретная модель при этом строится как первичная модель исследования, а способ перехода от континуума к системе с конечным числом степеней свободы основан, прежде всего, на физических соображениях.
Ниже описан численный метод анализа полей деформаций и напряжений, в котором в качестве дискретных моделей, представляющих сплошное тело, используются ориентированные графы [3-6]. Исследование системы на основе графового подхода сводится к тому, что 1) среда рассекается на отдельные части, имеющие известное математическое описание (в рассматриваемом случае — закон Гука); 2) для каждой части строится подграф (элементарная ячейка, графовый элемент), являющийся ее моделью; 3) элементарные ячейки объединяют в граф — модель анализируемого тела, после чего
с помощью матриц, характеризующих структуру графа, и уравнений, описывающих элементарные ячейки, получают уравнения системы в целом.
Графовая модель упругой среды в осесимметричной постановке рассматривалась в статье [3], где деформации в пределах элемента предполагались постоянными, и в статье [4], где используется линейная аппроксимация деформаций, однако при определении коэффициентов было сделано несколько упрощающих допущений. В данной работе при выводе матрицы жесткости осесимметричного графового элемента никаких упрощающих положений нет.
Основой для построения графов, используемых в качестве моделей, служит объективный характер операций измерения выбранных независимых переменных. Целью любого измерения является установление зависимости между значением переменной и показаниями прибора. В то же время измерение связано с точками системы, между которыми оно осуществляется, и ориентацией измерителя. Эти свойства операции измерения можно представить направленным отрезком, соединяющим две точки системы, т. е. дугой графа.
Таким образом, дискретную модель тела в виде ориентированного графа можно интерпретировать как стилизованную сеть приборов, установленных на тело для измерения его деформированного состояния.
При построении графа в качестве исходных используем параллельные переменные, т. е. переменные, которые измеряются установкой прибора непосредственно на тело [7]. Используя соответствующий прибор, например линейку, тензодатчик, тензометр и т.д., можно осуществить измерение линейных и угловых разностей перемещений, деформаций, положений точек относительно фиксированной системы координат.
Способ конструирования графа тела связан с процессом измерения полного и независимого комплекта параллельных переменных, которые однозначно характеризуют деформированное состояние элементов, полученных в результате декомпозиции. Элементарной ячейкой будем называть подграф, соответствующий одному элементу, полученному при разбиении исходной области на отдельные мелкие части.
При построении конфигурации элементарной ячейки в осесимметричной постановке сечение тела плоскостью, проходящей через его ось, покрываем сеткой координатных линий r _ const, z = const, между узлами которых устанавливаем гипотетические измерители. Полученную сеть измерителей преобразовываем в граф, представляя каждый из приборов дугой, а узлы, между которыми осуществлялось измерение, — вершинами графа. При этом вместо собственно измерителя будем ставить символ, отражающий природу измеряемой величины. Если вершины элемента обозначить A, B, C и D, то измерители необходимо установить так, чтобы измерять относительные перемещения точек A, B, C, D друг относительно друга (деформации 5rr, 5rz, 5zz, Szr) и относительно глобальной системы координат (перемещения
Ur , Ur , Ur , u D). Поскольку одни и те же точки участвуют в разных группах измерений — относительно осей r и z, это учитывается при построении графа тем, что соответствующие точки представлены разными вершинами. Ориентацию дуг, как и ориентацию процесса измерения, задаем в соответствии с направлением осей координат.
Каждой параллельной переменной, использованной для конструирования графа, соответствует парная ей последовательная переменная [7]. Выбор связной пары переменных определяется тем, что произведение последовательной и параллельной переменных должно давать скаляр с размерностью мощности или работы. Поскольку при построении графа в качестве параллельных переменных выбраны разности перемещений узловых точек либо перемещения узлов в глобальной системе координат, то каждой дуге графа, отображающей одновременно связную пару, соответствует последовательная переменная с размерностью силы. Точки приложения сил к элементу определяются концевыми вершинами дуги, т. е. узлами элемента.
В соответствии с описанной последовательностью получаем элементарную ячейку прямоугольного элемента в осесимметричной постановке, которая состоит из двух компонент (рис. 1). Дуги элементарной ячейки вводят два вектора — вектор деформаций {5} и вектор внутренних сил {f}:
_ {5d 5I u-D for uA u-B fod 5I 5u 5r }
ШТ _ r frr fzz fr fr frr fzz fr fr fzr frz fzr .erz\
e _ {fd ,fl ,fC, fD, fu ,fr ,fA, fB, fd ,fl , f u ,fr },
e
Т
e
(1) (2)
где S*r, ¿L, Кг, мации, , Slrz, S'Ur,
8rzz — нормальные дефор-относительные
перемещения, обусловленные поворотом сто-
,a
B
C
,D
Т
0OG
рон элемента, пГ, щ , иг , щ — перемещения точек, соответствующих вершинам в глобальной системе координат, /Гг, , /Гг, — нормальные внутренние силы (обобщенные напряжения), /|г, /р, , /Г2 — тангенциальные внутренние силы, /А, /В, /с, /Ь — внутренние силы в глобальной системе координат.
Индексы (п,й, 1,г) обозначают верхние (и) , нижние (й), левые (I) и правые (г) дуги элементарной ячейки, а индекс е указывает на принадлежность переменных отдельному графовому элементу.
Описание элементарной ячейки, т. е. связь векторов {/}е и {<5}е , а также зависимость их от напряжений и деформаций моделируемого тела устанавливаем, принимая в качестве инварианта перехода к модели энергию деформации элемента среды:
Рис. 1. Элементарная ячейка в осесимметричной постановке: а — компонента иг, б — компонента иг. Стрелками показано направление обхода контура, в кружках — верхний и нижний индексы переменных, представленных дугами ячейки
{f }T{¿}e = /МТ{в} dv,
(3)
где {<5}е и {/}е определяются формулами (1), (2). В результате энергию элементарной ячейки можно представить, с одной стороны, в виде
{f }T {^}e — fdr à?r + fZZ ^Zz + fC ur' + fD uD + fUr ^rr + frz +
i f r A + fr u B + fzr rd | frz ri | fzr ru + frz rr +fAur + JB ur + Jd °zr + Jl °rz + Ju °zr + Jr °rz,
(4)
с другой стороны, энергия деформации элемента сплошной среды определяется формулой
{a'}T{e} dv — (a/rr£rr + a/zzezz + a'**e^ + a/rzYrz) dv.
(5)
Используя дифференциальные зависимости Коши, выразим ковариантные компоненты тензора деформаций через перемещения
dur dr
duz du^
ezz ~ , eww ~ + rur 7
d^
dz
1 ( dur du
2 v "d^
1 f dur duz
erz ezr 7> I "Т; +
2 V dz dr
— e^r — 2 ( + dr
l / duz du^ — ^ — 2 ^^ + —
В случае осесимметричной задачи п^ =0 компоненты пг и п2 не зависят от а потому е^ = гпг, = 0, = 0. Закон Гука для изотропного тела запишем, используя тензор плотности напряжений {а'}, компоненты которого имеют вид а'7 = ^/дст- = га-. В результате получим
a/rr — (А + 2^)rerr + Arezz + , a/zz — Arerr + (A + 2^)rezz + ,
a'^ — Ae + Ae + A + 2Де
U C-rr I czz I Q
r r r^
/rz
a — ^r7rz,
(6)
где Yrz — erz + ezr, A и д — упругие постоянные Ламе, причем А — Ev/[(1 + v)(1 — 2v)] — для случая плоской деформации, А — Ev/(1 — v2) — для плоского напряженного состояния, д — 0, 5E/(1 + v), E — модуль упругости, v — коэффициент Пуассона.
r
Если ввести векторы напряжений [<r}T = {a'rr,a'zz ,a'^^,a'rz} и деформаций {e}T = {err, ezz, £^^,lrz}, то уравнения (6) можно записать в матричном виде:
{<} = [C]{e}, (7)
где [C] — матрица коэффициентов упругости.
Уравнения равновесия, записанные также с использованием тензорных плотностей, имеют вид [8]
da'rr da'rz d<'zz da 'zr
+ ^ - =0, ^ + ^ = 0. (8) dr dz dz dr
Для определения энергии деформации (5) необходимо проинтегрировать выражение
Irr- I _/zz _ I I _/rz.,
a £rr + a £zz + + a Yrz,
содержащее неизвестные деформации и напряжения. В связи с этим аппроксимируем неизвестные деформации в пределах элемента полиномами, линейными относительно переменных r и z:
dur duz т , , dur
£rr = = ao + ai r + a2 z, £zz = = bo + bi r + b2 z, = co + cir + C2 z,
dr dz dz
duz dur duz
-г— = do + dir + d2 z, Yrz = --+ -7т~ = Co + do + (ci + di )r + (c2 + d2)z, e= rur. (9)
dr dz dr
Заметим, что аппроксимирование деформаций, а не перемещений как в МКЭ, значительно сокращает размер матриц, используемых в дальнейших расчетах, а также исключает ряд проблем, возникающих при смещениях элементов как твердого тела. Выбирая локальную систему координат так, чтобы ось z совпадала с осью симметрии тела, а ось r проходила через центр элемента, в пределах элемента имеем Ri < r < Re, -0.5Az < z < 0.5Az.
Выразим теперь неизвестные коэффициенты полиномов в (9) через деформации сторон элемента. С помощью интегрирования вдоль соответствующей стороны элемента получаем
Re Az/2
/du f du
dr = (ao + aiRc + a2z)Ar, Srz = ~dz dz = (co + Cir)Az,
Ri -Az/2
Az/2 Re
/du f du
dz = (bo + bir)Az, özr = dr = (do + diRc + d2z)Ar,
-Az/2 Ri
где Ar = Re - Ri, Rc = 0.5(Ri + Re). В результате имеем:
S'Ur = (ao + aiRc - 0.5a2Az) Ar, 5f.r = (ao + aiRc + 0.5Az) Ar, Slrz = (co + ciRi) Az, Srz = (co + ci Re) Az, S'zz = (bo + bi Ri) Az, Srzz = (bo + biRe) Az, 8azr = (do + diRc - 0.5d2 Az) Ar, 5dzr = (do + diRc + 0.5d2 Az) Ar,
откуда получаем
Su + öd Sl + Sr öd - Su Sd- Su
rr rr rr - rr r - r
---, bo + bi Rc = ---, a2 = —:-:-, Ü2 = —:-:-,
2Ar ' o i c 2Az ' 2 ArAz ' 2 ArAz '
ör -Sl Sr -Sl Sl + ör su + sd
а = Srz Srz , bi = Szz Szz, co + ci Rc = rz + rz, do + di Rc = zrn + zr. (10) i ArAz ' i ArAz ' o i c 2Az ' o i c 2Ar v ;
Заметим, что
<*2 = С1, Ъх = d2. (11)
Это следует из контурного закона Кирхгофа [9], примененных к иг и и2 компонентам элементарной ячейки:
^ + & - ^ - ¿и = о, С + € - - К* = о
и соответствующих выражений из (10). Такой же результат следует из равенств
д2иг д2иг д2и2 д2и2 дгдг дгдг' дгдг дгдг'
В (10) найдена только часть неизвестных коэффициентов из (9). Чтобы найти оставшиеся, воспользуемся формулами (6), (9) и, обеспечивая повышенную точность метода, удовлетворим в пределах элемента уравнениям равновесия (8). В результате получим
г
Лй2Г + (Л + 2д) 62Г + Л (Со + С1Г + 02^) + д [со + ¿о + 2 (С1 + ¿1) г + (С2 + ¿2) г] = 0. (13)
(Л + 2д) (^ао + 2а1Г + а2г - ^Г^) + дг (с^ + ¿2) + Л61Г = 0, (12)
2д) 6
Из уравнения (12) следует
иг = аог + 2а1 г2 + а2гг + 1 [д (02 + ¿2) + Л&1 ] г2. (14)
Л + 2д
Кроме того, из (9) имеем
иг = ао г + 0.5а1 г2 + а2гг + / (г), (15)
где /(г) — произвольная функция г.
Из равенства выражений (14) и (15) для иг получим
/(г) = 1.5а1г2 + —Ц- [д(с2 + ¿2) + Л61 ] г2, (16)
Л + 2д
а это возможно, если 1.5(Л + 2д)а1 + д(с2 + ¿2) + ЛЬ1 = 0. Таким образом,
иг = аог + 0.5а1 г2 + а2гг. (17)
Уравнение (13) будет выполнено при условиях
Лсо + д(со + ¿о) = 0,
ЛС2 + д(с2 + ¿2) = 0, (18)
Ла2 + (Л + 2д)&2 + Лс1 + 2д(с1 + ¿1) = 0. (19)
Из (11), (16) и (18) находим, что
с2 =--^2 =--, (20)
д + Л д + Л
2Л
а1 = - 3(ЛТд)б1. (21)
В работе [4] при анализе уравнения (12) член а2г отбрасывается и принято дополнительное условие иг = (ао + а1Де)г , что не ведет к точному выполнению (12) в пределах элемента. Найдем оставшиеся коэффициенты 62 и . Из уравнения (19)
Ла2 + (Л + 2д)62 + (Л + 2д)с1 + 2д^ = 0 (22)
с учетом (11) получаем
Лс1 + (Л + 2д)62 + (Л + 2д)а2 + 2д^1 = 0. (23)
Сложив (22) и (23), находим
2(Л + д)а2 + 2(Л + 2д)62 + 2(Л + д)с1 + 4д^ = 0. (24)
Группируя в (24) коэффициенты, связанные с нормальными и тангенциальными составляющими, и приравнивая каждую группу к нулю, получим
Л + д Л + д
62 = - Л . 0 а2, ¿1 =---— с1. (25)
Л + 2д 2д
Используя (17), найдем перемещения вершин элемента ABCD: u^ = a0Ri + 0.5ai R2 — 0.5a2Ri Az,
uf = aoRe + 0.5a1R2 — 0.5a2ReAz, u^ = a0 Ri + 0.5a1R2 + 0.5a2 RiAz, u^ = a0 Re + 0.5a1 R2 +0.5a2 Re Az,
откуда получим
a0 + 0.5a1 Rc +
Ar2 4R
ur I u^- I 'uс I ur
4R
Таким образом,
u r = r
uA + uB + uC + uD 1 / „ Ar2 ,
-4R-+ 2a1 Г — Rc— 4R + a2z
Вводя в (9) для сокращения записи обозначения a = a0 + a1Rc, b = b0 + b1 Rc, c = c0 + c1 Rc, d = d0 + d1Rc, имеем:
err = a + a1 (r — Rc) + a2z, = b + b1 (r — Rc) + b2z,
Yrz = c + d + (c1 + d1 )(r — Rc) + (c2 + d2)z,
Ar2
(26)
ur I ur I ur I ur
4R
1
+ 2 a1
r — Rc — 4RT1 + a2 Z
где коэффициенты представлены в (10), (20), (21) и (25). В результате деформациям (26) можно придать вид
^ 2А
S^. + S£r 2Ar
3(А + д) ArAz
Sr - S1 Sd - Su
uzz "zz (r _ R ) + rr "rr
ArAz
z,
-
S' + 5r 5r - S'
2Az
+
ArAz
-(r — Rc) —
А + д — S,
z,
7rz —
S1 + sr Su + Sd
2Az
+
2Ar
+
д
— А Sr — S'
А + 2д ArAz
А Si - S
— r
ur I ur I ur I ur
2д А
ArAz
Sr - S1
zz — zz
(r — Rc) +
4Rc
3(А + д) ArAz
r Rc
А + д ArAz
Ar2 \
z,
(27)
4Rc
S^ — ¿ur + ArAz z
Уравнения (27) позволяют представить относительные деформации внутри прямоугольного элемента через абсолютные деформации его сторон
{е} = [L]{S},
(28)
где
[L] =
11 ¿12 0 0 ¿15 ¿16 0 0 0 0 0 0
21 ¿22 0 0 ¿25 ¿26 0 0 0 0 0 0
31 ¿32 ¿33 3 ¿35 6 3 ¿37 ¿38 0 0 0 0
0 0 0 0 0 0 0 0 ¿49 ¿4,10 ¿4,11 ,4 2
Элементы матрицы [L] имеют такой вид
¿11 = 0.5/Ar + z/(Ar Az), I12 = —¿16 = 2kb (r — Rc)/(3Ar Az), ¿15 = 0.5/Ar — z/(Ar Az), ¿21 = —¿25 = — kaz/(ArAz), ¿22 = 0.5/Az — (r — Rc)/(ArAz), ¿26 = 0.5/Az + (r — Rc)/(ArAz), ¿31 = —¿35 = r2 z/(ArAz), ¿32 = — ¿36 = kb r2 (r — Rc — 0.25Ar2/Rc)/(ArAz), I33 = ¿34 = ¿37 = ¿38 = 0.25r2/Rc, ¿49 = 0.5/Ar + k6 z/(ArAz), ¿4,11 = 0.5/Ar — k6 z/(ArAz), ¿4,10 = 0.5/Az — kc (r — Rc )/(ArAz), ¿4,12 = 0.5/Az + kc (r — Rc )/(ArAz),
где ka = (А + д)/(А + 2д), k6 = А/(А + д), kc = 0.5(д — А)/д. Подставив {е} из (28) в уравнение (3), получаем
{f}T {S}e = /{a'}T [L]{S}e dv = /{a'}T [L] dv {S}
е^^ - rur - r
¿
А А Тырымов. Графовый подход при построении конечно-элементной модели упругих тел Сократив последнее равенство на {5}е и осуществляя операцию транспонирования, имеем
а учитывая (7), получаем
{/}е = уМТ{*'} ¿V,
V е
{/}е = [^ [ОД{5}е¿V.
В результате уравнение состояния элементарной ячейки приобретает вид
{/}е = [К ]е {5}е,
где
[К]е = /МТ [ОД ¿V
V е
(29)
(30)
— матрица жесткости элементарной ячейки, элементы которой находятся в результате матричного умножения и последующего интегрирования. Уравнение (29) удобно представить в виде
/п
Л
Кп 0
0
К;
5п
где связь между нормальными составляющими {/п}, {5п}, [Кп] имеет вид
УаГ «11 «12 «13 «14 «15 «16 «17 «18
/ •1 22 «22 «23 «24 «25 «26 «27 «28
/с «33 «34 «35 «36 «37 «38
/ь «44 «45 «46 «47 «48
/и «55 «56 «57 «58
Л Г «66 «67 «68
/а «77 «78
_ /в . «88
;22
5иг 5Г
а связь между тангенциальными составляющими {/;}, {5;}, [К]; определяется следующим образом:
611
с »гг ^ /а
-/Т2 /г
/ 2Г /и
/ Г2 /Г
^12 &22
613
&23
Ьзз
&14
&24
Ьз4
644
Г 5а 2Г
5г Г2
< 5и Г
5Г Г
Элементы симметричных матриц [Кп] и [К;] таковы
«и = «55 = е1 Ле Аг/Аг, а12 = «25 = 0.25АЛе + е2 Аг, «13 = «14 = а17 = а18 = «35 = «45 = «57 = «58 = 0.125ААг, а15 = езДс Аг/Аг, «16 = «56 = 0.25АЯе - е2Аг, «22 = е4ЯсАг/Аг - (А + 2д)Аг2/(12Аг) + е5Аг3/(ЛсАг), «23 = «24 = «27 = «28 = 0.125ААг — е6Аг2/Ле, «26 = е7ЛсАг/Аг — е5Аг3/(Лс Аг),
«33 = «34 = «37 = «38 = «44 = «47 = «48 = «77 = «78 = «88 = (А + 2^)АгАг/(16#с),
«36 = «46 = «67 = «68 = 0.125ААг + е6 Аг2/Ле, «66 = е4 Яс Аг/Аг + (А + 2д)Аг2/(12Аг) + е5 Аг3/(ЯсАг), Ьц = 633 = 0.25^(1 + е8)Дс Аг/Аг, 612 = 623 = 0.25д(Яс — ед), 613 = 0.25^(1 — е8 )ЯС Аг/Аг, 614 = 634 = 0.25^(ЯС + ед), 622 = 0.25д[Яс(1 + ею) — 2ед ]Аг/Аг,
е
е
Ь24 = 0.25^(1 - е10)ЛСДг/Дг, Ь44 = 0.25д[Де(1 + е10) + 2е9]Дг/Дг,
где
= 4А2 + 22Лд + 21д2 = Л(д - 3Л) = 2Л2 + 2Лд + 3д2
е1 = 12(Л + 2д) ' е2 = 72(Л + д)' ез = 12(Л + 2д) ' = 27Л3 + 136Л2 д + 180Лд2 + 72д3 = Л2 (3Л + 8д) = Л(3Л + 7д)
е4 = 108(Л + д)2 ' е5 = 432(Л + д)2 ' е6 = 144(Л + д)'
= 27Л3 + 80Л2д + 90Лд2 + 36д3 = Л2 = д - Л Д = (д - Л)2
е7 = 108(Л + д)2 ' е8 = 3(Л + д)2' е9 = 12д"' е10 = 12д2 '
Таким образом, для четырехузлового элемента получена матрица жесткости при аппроксимации поля деформаций линейными полиномами. Стандартный МКЭ для этого требует применения элементов с восемью узлами. В результате при таком же уровне аппроксимации графовый метод использует систему уравнений, содержащую примерно в 3 раза меньше уравнений по сравнению с системой, полученной традиционным способом в МКЭ. Тело, представленное в виде отдельных элементов, и соответствующая ему совокупность элементарных ячеек описывается уравнением
{/} = [КМ, (31)
где{/} и {5} — векторы внутренних сил и деформаций
{/}Т = {/Те, /Те, ■■■, /Те}, {5}Т = {^е, , ■■■, },
причем {/}^е, {5}^ определяются (1), (2), а [К] — глобальная несвязанная матрица жесткостей тела, представленного в дискретном виде
[К] [К1е,К2е, ■ ■ ■ КПе],
где [Копределены в (30), а п — число ячеек, образующих граф.
Уравнение (31) связывает внутренние силы с деформациями элементов. В общем случае с помощью этого уравнения деформированное состояние определить нельзя, поскольку, как правило. заданными бывают внешние силы и перемещения. Теория графов обеспечивает возможность искусственного конструирования квадратных матриц преобразования, с помощью которых можно, исходя из уравнений, описывающих разрезанное на элементы тело, получить уравнение связного тела. Структурные элементы графа — разрезы, контуры, пути и хорды — могут быть заданы матрицами. Примеры этих матриц для элементарной ячейки (рис. 1) приведены в [4] .
После преобразования с помощью квадратных матриц дополнений [Ави маршрутов [Рсаналогично [4,6] имеем
{Е Ь = [Ав Ь {/}, {ДЬ = [Рс Ь {5}, (32)
где индексом Б обозначены переменные глобальной системы координат графа.
Поскольку матрицы дополнений и маршрутов несингулярны, а также в силу равенств [Рс]Т = [Ав]-1 и [Ав]Т = [Рс]-1, получим
{/} = [Ав ]-1{^ Ь = [Рс ]£ {^ Ь, {5} = [Рс ]-ЧДЬ = [Ав ]£ {ДЬ ■ (33)
Тогда, используя (31), (32), (33), получаем уравнения, описывающие связное тело
{Е ь = [К^ {ДЬ, (34)
где [К]з = [Ав[К][Ав— глобальная матрица жесткости тела.
Использованный для вывода уравнения (34) подход, можно рассматривать как модифицированный метод конгруэнтных преобразований [2].
Для иллюстрации работы метода решим несколько тестовых задач.
Задача 1. Рассмотрим задачу изгиба круглой однородной пластины радиуса Л, защемленной по контуру. Пластина нагружена равномерно распределенным внешним давлением Р/Е = 0-001. Сравнение результатов расчетов проводилось путем сопоставления решений, полученных на достаточно
редких сетках с использованием различных конечных элементов [10], элементарных ячеек графового метода, а также точного решения. В табл. 1 приведены значения прогиба в центре пластины толщины h/R = 0.1, выраженные в долях PR4D-1, где D = Eh3/(12(1 — v2)) — цилиндрическая жесткость пластины. Значение коэффициента Пуассона принято равным 0.3. Использовались треугольные конечные элементы Морли с шестью степенями свободы [11], BCIZ (элемент Бейзли, Ченга, Айронса и Зенкевича) с девятью степенями свободы [12], элемент DKT (discrete Kirchhoff triangle) с девятью степенями свободы [13], гибридный метод напряжений HSM с девятью степенями свободы [10]. Эти треугольные элементы построены для тонких пластин, рассчитываемых по теории Кирхгофа. В расчетах графовым методом применялся как четырехугольный элемент, представленный в данной статье, так и треугольный элемент графовой модели в осесимметричной постановке [5].
Таблица 1
Элемент Элемент Элемент Элемент Треугольный Четырех- Точное Точное
Морли BCIZ DTK HSM элемент угольный решение решение
графового элемент по теории по теории
метода графового Кирхгофа Тимошенко-
метода -Миндлина
0.01914 0.01643 0.01597 0.01541 0.01565 0.01628 0.01563 0.01634
Во всех случаях количество элементов, на которые разбивалась область, принято равным 96.
Задача 2. Рассматривается осесимметричная задача о напряженном состоянии в полом цилиндре под действием нормального давления, равномерно распределенного по внутренней боковой поверхности (задача Ламе). Аналитическое решение для вектора перемещения вдоль радиуса имеет вид
иг = г?(1 + V)Е-1 р[(1 - 2v)г + г?/г]/(г? - г?),
где г1, г? — внутренний и внешний радиусы цилиндра соответственно, р — равномерное внутреннее давление.
В расчетах принято: г1 = 3 мм, г? = 7мм, р = 1 МПа, Е = 100 МПа, V = 0.3. На торцах в выделяемом для анализа фрагменте цилиндра требовалось условие иг = 0. В табл. 2 приведено сравнение теоретических и расчетных результатов для перемещения иг в зависимости от радиуса при различных сетках, покрывающих рассматриваемую область. Во второй и третьей строках табл. 2 даны расчетные значения, а также их процентное отношение к теоретическим для сетки 9 х 6. Результаты в четвертой и пятой строках относятся к сеткам 90 х 6 и 90 х 60 соответственно.
Задача 3. Рассмотрим осесимметричную задачу о напряженном состоянии сплошного изотропного цилиндра высотой 2Н со свободной от напряжений боковой гофрированной поверхностью, который находится под действием осевых растягивающих усилий интенсивностью р (рис. 2). В безразмерных цилиндрических координатах, отнесенных к внешнему радиусу невозмущенного цилиндра, гофрированная поверхность описывается уравнением г = 1 + е соб пг.
Ввиду симметрии рассматривается четверть всей расчетной области. На осях симметрии должны выполняться следующие условия: иг = 0 при г = 0, иг =0 при г = 0 (при соответствующих значениях другой переменной).
Числовые расчеты проведены для изотропного цилиндра высотой 2Н = 4 с коэффициентом Пуассона V = 0.33 при е = 0.1. На рис. 3 показано распределение напряжений по толщине цилиндра в сечении г = 1. Напряжениям агг/р, /р, и99/р соответствуют пунктирная, сплошная и штриховая линии.
Таблица 2
r, мм 3 4 5 6 7
ur, мкм 51.29 40.51 34.52 30.91 28.67
Ur /игтеор 100.000 99.8836 99.4679 99.9432 100.000
% 100.000 99.9947 99.9976 99.9978 100.000
% 100.000 99.9948 99.9978 99.9978 100.000
■ОТ V V V V WW
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Рис. 2. Расчетная схема
Рис. 3. Распределение напряжений при различных значениях г в сечении г = 1
Сравнение этих результатов с решением данной задачи полученным методом возмущения формы границы [14] дает хорошее совпадение.
Тестирование показывает весьма высокую точность результатов, получаемых на основе предлагаемого метода. Повышенная точность вычислений обеспечивается благодаря тому, что, с одной стороны, графовые закономерности реализуют выполнение условий равновесия и неразрывности деформаций для элемента в целом, а с другой стороны, уравнения равновесия выполняются локально по объему элемента.
Библиографический список
1. Флетчер К. Численные методы на основе метода Га-леркина. М. : Мир, 1988. 352 с.
2. Галлагер Р. Метод конечных элементов. Основы. М. : Мир, 1984. 428 с.
3. Кузовков Е. Г., Тырымов А. А. Графовая модель упругой среды в осесимметричной постановке // Моделирование в механике : сб. науч. тр. Новосибирск : Изд-во СО АН СССР, 1990. Т. 4(21), № 6. С. 103-109.
4. Kuzovkov E. G. Axisymmetric Graph Model of an Elastic Solid // Strength of Materials. 1996. Vol. 28, № 6. P. 470-485.
5. Тырымов А. А. Треугольный элемент графовой модели для осесимметричной задачи теории упругости // Численные методы решения задач теории упругости и пластичности : тр. XVIII Межресп. конф., Кемерово, 13 июля 2003 г. / под ред. В. М. Фомина. Новосибирск : Нонпарель, 2003. С. 187-191.
6. Тырымов А. А. Сингулярный элемент графовой модели упругой среды в декартовой системе координат // Вычислительная механика сплошных сред. 2011. Т. 4, № 4. С. 125-136.
7. Trent H. Isomorphism between Oriented Linear Graphs and Lumped Physical Systems // J. of the Acoustical
Soc. of America. 1955. Vol. 27, № 3. P. 500-527.
8. Крон Г. Исследование сложных систем по частям
— диакоптика. М. : Наука, 1972. 542 с.
9. Свами М., Тхуласираман К. Графы, сети и алгоритмы. М. : Мир, 1984. 454 с.
10. Белкин А. Е., Гаврюшин С. С. Расчет пластин методом конечных элементов. М. : Изд-во МГТУ им. Н. Э. Баумана, 2008. 232 с.
11. Morley L. S. D. The constant-moment plate-bending element // J. Strain Anal. 1971. Vol. 6, № 1. P. 20-24.
12. Bazeley G. P., Cheung Y. K, Irons B. M, Zienkiewicz O. C. Triangular elements in plate bending
— conforming and non-conforming solutions // Proc. Conf. on Matrix Methods in Structural Mechanics. Ohio : Wright-Patterson Air Force Base, 1965. P. 547-576.
13. Batoz J. L, Bathe K. J., Ho L. W. A study of three-node triangular plate bending elements // Intern. J. for Numerical Methods in Engineering. 1980. Vol. 15. P. 1771-1812.
14. Немиш Ю. Н. Элементы механики кусочно-однородных тел с неканоническими поверхностями раздела. Киев : Наук. думка, 1989. 312 с.
