книги / Напряженное состояние и прочность оболочек из хрупких неметаллических материалов
..pdfБудем предполагать, что тело является цилиндрически ортотропным. В этом случае через каждую его точку проходит плоскость, все направления которой эквивалентны в отношении упругих свойств.
Меридиональное сечение R тела вращения разобьем на ряд контак тирующих между собой треугольных элементов. Вершины этих треу гольников, обозначенные для каждого элемента /, /, k, являются узла ми области. Тело вращения, таким образом, разбивается на ряд торо образных элементов с сечением в виде треугольника (рис. 1, а).
Рис. 1. Разбиение области на треугольные элементы:
а — треугольный элемент lt /, k в осевом сечении тела вращения; 6 — задание поверхностных нагрузок в угловой точке i\ о — «звезда» из треугольных эле ментов вокруг узла /, в которых анализируются значения углов (J; г — схема задания поверхностных нагрузок в узле /.
Задача, записанная относительно непрерывных переменных (пере мещений), заменяется задачей относительно дискретных значений этих перемещений в узлах треугольников, Ur = Ut\ Uz = V{. (В дальней шем будет использоваться матричная форма записи для основных соот ношений МКЭ, что делает изложение и запись более компактными.)
Далее ход рассуждений будет вестись применительно |
к треуголь |
ному элементу i, /, k (рис. 1). |
элемента R\ |
Введем вектор перемещений в узлах треугольного |
|
; |
(П-1) |
Примем, что компоненты вектора перемещений внутри элемента распределяются линейно в зависимости от узловых перемещений
= |
/ = alf + <x2fz + о у , |
(II.2) |
где « 1/, <*2/, аз/ могут быть найдены из двух систем совместных уравне ний, которые получаются из каждого уравнения (II.2), если вместо
г, г подставить координаты узлов, а перемещения f |
приравнять уз |
ловым: |
|
/ = "2д- [(at + М -f с(г) ft -f (й/ + bjZ -{- c,r)fi + |
|
+ (ak + bkz + ckr)fk). |
(II.3) |
Здесь A — площадь треугольника i, j, k\ |
|
at = z f k— zkrj\' |
|
b i ^ r , — rk\ |
(II.4) |
cl = zk —
Другие коэффициенты получаются в результате циклической пере становки индексов. В матричной форме соотношения (II.4) при мут вид
|
|
|
(II.5) |
где I — единичная матрица (2 X 2); |
|
|
|
П‘ ~ |
2Д |
* |
(П.6) |
IN], IN], INI — матрицы функций |
«положения», которые дают соот |
||
ветствующие узловые перемещения, если в уравнение (II.5) |
подставля |
ются координаты соответствующих узлов.
Изменение перемещений в пределах элемента по линейному за кону (П.З) удовлетворяет условиям непрерывности в пределах всего тела. Вдоль каждой стороны треугольника перемещения изменяются линейно и, следовательно, полностью определяются значениями в двух соответствующих узлах. Поэтому для соседних элементов с одинако выми перемещениями в узлах перемещения в любой точке их общей стороны совпадают, и непрерывность перемещений не нарушается при переходе от элемента к элементу.
При выборе начала координат в центре тяжести треугольного эле
мента выражения для ah at, ak (II.4) изменяются! |
|
||||
at = |
а, = |
аь = 2Д/3; |
(П.7) |
||
|
|
1 |
г, |
|
|
2А = |
Det |
1 |
г, |
г, |
(П.8) |
|
|
1 |
г* |
гк |
|
Для осесимметричной задачи вектор деформации, определенный че рез перемещения узлов, можно представить в виде
|
дг |
|
|
®г |
dW |
|
|
®2 |
дг |
|
(П.9) |
U |
|
||
6ф |
|
||
|
|
||
, &kz, |
|
|
|
bU |
. |
dW |
|
дг |
"** |
дг |
|
Деформации можно представить в матричной записи |
|||
{е}я = |в]. м * = [В,в а ] • |
(Н.ю) |
||
где |
|
|
|
°1 |
|
0 - |
|
0 |
|
bt |
|
в{ = |
|
0 |
; |
•у + “ + |
|
||
*i |
|
0{ _ |
|
B/t Вk можно получить в результате циклической перестановки ин дексов.
Матрица [Б] включает координаты г, 2, является их функцией, поэтому деформации в пределах элемента не будут постоянными. Иногда в теле возникают так называемые начальные деформации, выз ванные усадкой, ростом кристаллов, температурными изменениями, и не зависящие от напряжений. Обозначим их
БгО
Ы R |
6z0 |
(ii.li) |
£ф0 |
||
|
£/*zo |
|
Д ля изотропных материалов температурная деформация, вызванная неравномерным нагревом, будет представлена компонентами вектора
а Г 0
R |
«Т’о |
( 11. 12) |
{®о} |
а Г 0
Здесь Та — средняя температура элемента; а — коэффициент темпе ратурного расширения.
В предположении упругого характера работы системы связь между напряжениями и деформациями в виде закона Гука записывается в матричной форме
|
°г |
' |
|
|
{о} = |
Ог |
= [D] ({е} - {е0}), |
(11.13) |
|
Оф |
||||
|
|
|
>^гг ■
где [D] матрица упругости.
Построим матрицу упругости для цилиндрически ортотропного те ла. Известно, что в этом случае имеется пять независимых «технических» постоянных: Ег и £ 2— модули Юнга соответственно для направлений в плоскости изотропии (плоскости Гф) и для направлений, перпенди кулярных к ней; Vj и v2 — коэффициенты Пуассона, соответственно характеризующие сокращения в плоскости изотропии при растяжении
в направлении, перпендикулярном к ней; |
G2 — модуль сдвига, |
харак |
||||
теризующий искажение |
углов |
между направлениями в плоскости |
||||
изотропии и направлением, перпендикулярным к ней. |
|
|||||
Обозначим |
Ег |
|
|
|
|
|
|
= я; |
G« |
т; |
(11.14) |
||
|
У |
у |
= |
|||
тогда |
|
|
|
|
|
|
~ 1— V* nv2(l + |
vx) |
nv2(1 + |
vx) |
0 |
|
|
« ( 1— nv*) |
«Vj + |
Vi |
0 |
|
||
D = а |
|
n (1 — rtVp |
0 |
.(11.15) |
||
Sym |
|
|
|
|
m( 1 + vx) x |
|
|
|
|
|
|
X (1 — vx — 2/IV2) _ |
|
Здесь |
|
|
|
|
|
|
а = |
(1+ vt) <1—Vj — 2т ф |
|
Для изотропного материала
I |
I |
- |
Gj |
|
|
Е* |
= 1; |
|
|
и матрица (11.15) преобразуется к виду |
||||
|
|
|
1 |
V |
|
|
|
1—V |
|
D = |
|
E(1- v ) |
|
1 |
U + v ) ( l - 2v) |
|
|
Sym
2(1 |
1 |
(11.16) |
|
+ v) ’ |
|||
|
V
1- v 0
V
1—v 0
(11.17)
1 0
1—2v 2(1 — v)
Внешние нагрузки, объемные и поверхностные, а также заданное поле перемещений соответствует некоторым узловым статически им эквивалентным нагрузкам, которые можно определить, используя на чало виртуальных перемещений. При этом задаются произвольные
(виртуальные) |
перемещения и приравниваются работы внешних сил |
и статически |
им эквивалентных узловых нагрузок на этих виртуаль |
ных перемещениях. |
|
Введем вектор обобщенных сил (Р}я в узлах элемента R, отвечаю |
|
щий вектору объемных сил (0}Л, вектор обобщенных узловых сил { P } \ |
отвечающих заданным поверхностным силам {F}*, и вектор обобщен
ных узловых усилий (5}Л, статически эквивалентный заданному полю перемещений. Зададим в узлах виртуальное перемещение (б). Тогда приращение перемещений и деформаций внутри элемента соответствен но имеют вид
|
|
S i f f = |
IN] б {?}*; |
|
|
|
|
|
б{е}* = |
[В]б{<7}*. |
|
(II-18) |
|
Приравняв соответствующие работы |
векторов |
(0}Л и [ P f, |
{F}* и |
|||
(Р)Л на |
виртуальных перемещениях, |
получим |
|
|
||
|
I |
(M 6{< 7)V {0}*dK = |
J (6 {?}У |
|
(И.19> |
|
|
уR |
|
|
vR |
|
|
|
J |
т б (<7}Л)Т{F}^ dS = |
j (б [ q f f |
{Р}* dS. |
(Н.20> |
|
|
|
|
|
s 1* |
|
|
Здесь ( |
)т обозначает операцию транспонирования. |
|
Учитывая, что объемный и поверхностный интегралы берутся повсей длине окружности, получаем
{Р}* = |
2я J |
J [Nf {Qfrdrdz; |
(11.21) |
|
|
SЯ |
|
{Р}* = |
2я J |
J [N)T{F frd rd z . |
(11.22) |
Для нахождения вектора (5}Л, статически эквивалентного задан ному полю перемещений, приравняем приращение работы внутренних, усилий на виртуальных перемещениях, т. е. приращения потенциаль
ной энергии деформации, и работу сил {5 }я на возможных вариациях, перемещений в узлах:
J ({«»)*)’ {S}*dK- J ({6б)яГ (6}dv =
yR |
|
yR |
- |
№)*)’ |
{«)*<»'- П ( { & 7 ) * ) W - |
yR |
|
\y R |
- |D1 |BJ dV\ {q\n = 2n у (((в?)'*)’ |fil' [D] [fij rdrdl) {»)*. |
(11.23) |
Таким образом, |
|
|
|
{ S f = |
(2я J J |
[В]т [D] [В] rdrdz) {q)R. |
|
Квадратная матрица |
|
|
|
[K}* = |
2J I J J |
[B]T[D]\B] rdrdz |
(11.24) |
называется матрицей жесткости. Она является блочной матрицей, элементы которой в свою очередь представляют подматрицы
|
Г/С« |
Кц |
К » ' |
|
|
К , |
К» |
Кр, , |
(Н.25) |
|
_Кы |
Кщ |
Kkk. |
|
где Kmi — квадратная |
(2 X 2) матрица, |
|
|
|
УСт|] = |
2я J J [В,]т [D] [BJ rdrdz. |
(11.26) |
||
Уравнение (11.24) примет вид |
|
|
|
|
|
{S|* - уд* (,)*. |
(11.27) |
Элементы матрицы [В] являются функциями координат г, z (II.Щ , поэтому получение простого выражения для матрицы путем интегри рования (11.26) сопряжено с большими трудностями. Приближенный способ для упрощения матрицы жесткости заключается в замене зна чений всех функций, входящих в матрицу 1В], значениями их в центре
тяжести треугольника |
г — Г{ ^ Tk-\ z — Z ^ ~^2-. Тогда |
при |
ближенная формула |
для матрицы жесткости запишется в виде |
|
|
[/С]=2я[В]т[В>][В]7д, |
(11.28) |
где [В]— матрица [В], вычисленная для центра тяжести треугольника;
А— площадь треугольника.
Вформулах (11.21), (11.22), (11.24)
\РЛ*i
{P)R = |
P, |
(11.29) |
|
P k |
|
[P)R = |
Pi |
|
Pj |
(11.30) |
IP-,S= 0
p k.
[Si
{S}* - |
(11.31) |
В случае действия на тело вращения объемных сил выражение (11.28) можно упростить, приняв начало координат в центре тяжести треу гольника.
При этом
j rdrdz — § zdrdz.
Используя выражения (II.6) и (II.7) и подставляя их в уравнение (11.21), получаем видоизмененное уравнение (11.29):
(11.32)
р0г
( р ) я = ) л д/3-
Формула (11.32) означает, что все объемные силы распределены в узлах поровну.
Перейдем к выводу разрешающих уравнений равновесия в пере мещениях. Эти уравнения можно получить, исходя как из начала воз можных перемещений, так и из вариационного принципа Лагранжа о стационарности полной потенциальной энергии тел, которая на основании формулы Клайперона может быть представлена в виде
П = |
П — А, |
(11.33) |
где П — потенциальная энергия |
деформаций, |
А — работа внешних |
поверхностных и объемных сил. |
|
|
Используя представления (11.18) — (11.24), выражение (11.33) пере пишем в виде
- J S |
(Р}я rdrdz - \ ((«*)’ [ P f rdl |
|
S |
■ |
l |
= 2я (0,5 |
$ I (?V |
[B]T|D| [B] {qfrdrdz - |
|
s |
|
- J J((?)Y vn’ (P)R rdrdz) = $ (to)*)T|WJT[ P f rdl
S |
l |
- 0.5 <(*)Y (S)S - ((«)Y (PI* - ((»}Y |P)r. |
(11.34) |
Вводя вектор внешней нагрузки в |
узлах |
= {P}R + {P\R |
и учитывая уравнение (11.27), переписываем выражение (11.34) |
||
П = 0,5 ([д\У I K f { q f - |
({д}*)т {Р}* |
(11.35) |
Известно, что в положении равновесия потенциальная энергия стационарна. Введем [£ х] диагональную матрицу с числом элементов по диагонали, равным порядку вектора {q}, причем каждому элемен
ту матрицы, стоящему на диагонали, |
поставим |
в соответствие |
узел |
и направление перемещений в том же |
порядке, |
что и в векторе |
{<7}. |
Если перемещения в узле заданы, соответствующий диагональный член принимается равным нулю, в противном случае диагональный член равен единице.
Условие стационарности потенциальной энергии на возможных вариациях вектора перемещения в узлах дает
41^ - °. <"-зв>
что в развернутом виде представляется следующим образом: |
|
IPiHIKl ~ {q} ~ {Р}) - 0. |
(11.37) |
Уравнение (11.37) является разрешающим уравнением равновесия
вперемещениях, полученным из условий стационарности и записанным
вматричной форме.
Раскроем смысл всех членов, входящих в уравнение (П.37): {qj — вектор перемещений всех N узлов тела; [£] — вектор внешней нагруз ки в узлах,
|
|
|
|
{ЯЛ |
) |
|
|
|
|
|
|
|
[R) |
{**} |
|
|
|
(11.38) |
|
|
|
|
|
|
|
|
|
||
|
|
|
i |
{Я *} |
|
|
|
|
|
{£<) — /-и блочный компонент, равный |
сумме векторов узловых |
сил |
|||||||
{P[)R по всем элементам, сходящимся в узле t, и действующих на |
эле |
||||||||
мент R со стороны узла i, |
|
|
REra& |
|
|
|
|||
|
|
{ R , ) = £ , < « ' > * = |
|
(11.39) |
|||||
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|||
Символ R € 1 означает |
суммирование по всем элементам, сходя |
||||||||
щимся в узле i; IК\ — матрица жесткости всей системы, |
|
|
|||||||
|
|
|
|
о |
|
Б |
l^ClAr]^ |
|
|
|
Б [Knf |
S |
[Kuf |
|
R £IN |
|
|
|
|
|
леи |
|
R£U |
|
|
|
Ww]R |
|
|
1/С1 = |
Б |
№ 1 * |
Б 1Kuf |
|
S |
(11.40) |
|||
|
R£IN |
|
|
|
|||||
|
R£li |
( K N I ]R |
R£U |
|
|
R£NN |
|
|
|
|
Б |
Б iW* |
|
|
|||||
|
|_-яем |
|
R£Nl |
|
S |
1К Ы * |
|
|
|
|
|
|
|
|
|
|
Здесь |
символ R £ iN — суммирование |
по элементам, содержащим |
узлы |
i и N. |
X 2) матрица, [KI N \ = 2я X |
Напомним, что [/С*л/] — квадратная (2 |
X ([Я/Г [D] 1Яуу])гД вычисленная для элемента R.
Подматрицы матрицы [/С;дг]— нулевые, если по крайней мере один из узлов i или N не принадлежит элементу R.
Решив уравнение (11.36), относительно компонентов вектора уз ловых перемещений {q}, можно по формулам (11.10) и (11.13) найти деформации и напряжения и таким образом полностью решить краевую задачу.
Так как матрица [В], как уже было сказано, зависит от координат, напряжения изменяются по всему элементу (11.13). Обычно вычисляют напряжения в центре тяжести треугольников, а затем усредняют их, суммируя для всех элементов, содержащих узел i.
Переход от функциональных уравнений теории упругости к их раз ностным аналогам выдвигает проблему решения систем линейных алгебраических уравнений высокого порядка. Методы решения спе цифических систем линейных алгебраических уравнений, порождае мых методом конечных элементов, постоянно совершенствуются [22, 54, 56, 87, 90, 98, 109].
В некоторых универсальных пакетах прикладных программ широ ко используется метод сопряженных градиентов, однако наиболее ши роко применяются различные модификации метода исключения Гаус са [75, 90, 107, 109, 145, 146].
За последние годы создан ряд эффективных программ, ориентиров ванных на решение плоских и осесимметричных задач теории упругос ти [56, 59, 60, 75, 109, 115]. Усилия разработчиков пакетов приклад ных программ сосредоточиваются на повышении степени автоматиза ции всех этапов решения задачи: от разбиения области заданного конструктивного элемента до выдачи графической информации о рассчи танных компонентах тензора напряжений и перемещениях, характе ризующих напряженно-деформированное состояние элемента от воздействия заданных нагрузок. Ниже кратко опишем методику ис пользования метода конечных элементов для решения упругих осесим метричных задач при помощи одного из высокоавтоматизированных пакетов прикладных программ «Элемент» [60].
Основная часть затрат инженерного труда и машинного времени при решении задач математической физики методом конечных элемен тов связана с подготовкой исходной информации и обработкой полу чаемых в результате расчета числовых массивов.
Предложенный алгоритм и соответствующая ему программа [60, 71 ] позволяют в значительной степени автоматизировать указанные выше операции. При этом производительность исследователя повыша ется в десятки раз, а также существенно снижается вероятность ошиб ки в исходных данных вследствие многократного уменьшения их объема.
Важнейшими условиями при разработке было обеспечение достаточ но высокого качества сетки независимо от специфических особенностей каждой конкретной области. К числу таких особенностей, затрудняю
щих автоматизацию, относятся большая кривизна; сильная вытянутость в одном направлении, наличие полостей, необходимость сгущения в ок> рестности концентраторов и т. д. При оценке качества построенной сетки следует принять во внимание, что форма и размеры элементов оказы вают существенное влияние на точность решения задачи, а затраты машинного времени определяются (при использовании прямых методов для решения систем линейных уравнений) топологической структурой сетки и способом нумерации узлов.
В связи G этим основные требования, предъявляемые к конечно элементному разбиению, следующие:
размер элементов в зонах ожидаемых больших градиентов искомой функции должен быть достаточно мал для того, чтобы обеспечить ее адекватную аппроксимацию;
нумерация узлов должна давать близкую к минимально возможной ширину ленты в системе линейных уравнений (ширина ленты пропор циональна максимуму разности номеров двух узлов, входящих в один и тот же элемент);
число узлов и элементов не должно превышать пределов, устанав ливаемых объемом памяти, доступной пользователю ЭВМ.
Удовлетворить сразу всем перечисленным выше требованиям для достаточно сложной области весьма трудно. В связи с этим в ряде случаев приходится прибегать к повторному построению сетки после коррекции исходных данных.
Перечислим основные этапы построения сетки при использовании автоматической триангуляции.
г-' 1. Разбиение контура области на отрезки и занесение в массивы координат г и г контурных узлов и номеров контурных узлов, следую щих в порядке обхода области против часовой стрелки.
2. Триангуляция области, т. е. подразделение ее на конечные тре угольные элементы.
3.Оптимизация нумерации узлов с целью снижения ширины ленты
вматрице системы линейных уравнений.
4.Пересортировка выходных информационных массивов, харак теризующих разбиение области с учетом изменившейся нумерации узлов.
Все этапы, начиная со второго, выполняются программой автомати чески. Первый этап представляет собой подготовку исходных данных к программе. При этом область может быть разбита на подобласти, границы которых включают линии раздела между различными материа лами [7, 24, 98]. Каждая из линий — границ подобластей — подраз деляется на отрезки, которые на следующем этапе работы алгоритма должны стать сторонами конечных элементов. Поэтому разбиение гра ницы в значительной степени предопределяет густоту и другие свойства сетки. В процессе подразделения линий координаты деления заносятся в массив координат узлов, а информация о граничных условиях, задан ная применительно к линиям, преобразуется и приводится к узловым точкам. Области и подобласти произвольной формы триангулируются последовательным заполнением области, предварительно подразделен ной на отрезки от Гранины во внутрь ее. Область заполняется так, что