Применение метода конечных элементов
..pdf12.1. Теория упругости. Одномерный случай
Простейшая одномерная задача является удобной отправной точкой для дальнейшего изложения, поскольку только она позво ляет проиллюстрировать все выкладки для конкретного числового примера. В двумерном и тем более в трехмерном случае объем вычислений слишком велик для этого. Так как обсуждение, прово димое ниже, относится к отдельному элементу, верхний индекс в обозначениях всех величин, за исключением [Л(<7)] и будет опускаться.
Предполагая, что одномерное упругое тело ориентировано вдоль оси х, будем иметь только одну компоненту тензора напря жений Охх и соответствующую компоненту тензора деформаций ехх. Запишем закон Гука
ахх—Е (Iе) — !ео))- |
(12.3) |
Формула (12.3) в матричном виде записывается как |
{о} = |
= [£ ]({« } — {е0} ), поэтому [D]=E, где Е — модуль упругости. На чальную деформацию обычно связывают с тепловым расширением аАТ, где а — коэффициент теплового расширения, а АТ— отклоне ние температуры от некоторого равновесного значения.
Для одномерного элемента функция перемещения имеет вид
« = е д + е д = [ л п ( ^ ) , |
(12.4) |
|||
где Ui и Uj — перемещения узлов |
i и / в направлении |
оси х. Д е |
||
формация ехх связана с перемещением формулой |
|
|||
du |
dNi |
dNj |
(12.5) |
|
&хх dx |
dx |
dx |
||
|
||||
Производные от функций формы вычисляются легко, так как |
||||
* . = 1 - Т |
» * /= - £ " |
|
||
Дифференцирование дает |
|
|
|
|
e « = - f l - l |
Ц \Щ = { В ][ и ). |
(12.6) |
Матрица градиентов [В] теперь определена, так что можно при ступить к составлению матрицы жесткости. Подставляя [В] и [D] в формулу (12.1) и предполагая площадь поперечного сечения по стоянной, получаем
|
L |
|
|
[кЩ = [В\т[D] [ B \ d V = ^ t |
— 1 |
dx, |
(12.7) |
[ - 1 I j f |
1
ИЛИ
т = т |
1 |
( 12. 8) |
|
|
— 1 |
Соотношение (12.8) идентично по форме матрице элемента, полу ченной в одномерном случае переноса тепла.
Интегралы, определяющие вектор нагрузки, вычисляются так же просто. Интеграл, связанный с тепловым расширением, записы вается как
- j W [D] {е0} d V = - Щ т |
|
J - jJ j* dx= - а Е А (AT) J“ |
JJ. (12.9) |
|
|
|
L |
|
|
v |
|
о |
|
|
Интеграл от объемных сил имеет вид |
|
|
||
- J w m d V = —XA |
оI |
dx= |
гI A L (i |
( 12. 10) |
|
||||
V |
|
|
|
Из поверхностных нагрузок в одномерном случае остается только рх, и она должна быть сосредоточена в одной из узловых точек. Предполагая, что рх приложена в t-м узле, вычислим поверхност ный интеграл
(12.11)
где нижний индекс / обозначает номер узла. Матрица функций формы сводится к { J }, потому что нагрузка сосредоточена в узле. Если она приложена в /-м узле, поверхностный интеграл записы вается как
- J Px[ jv r d s = - M , { J } - |
02.12) |
s
Полная система уравнений, определяющих элемент, имеет вид
(12.13)
Все объемные интегралы должны быть вычислены заново, если площадь поперечного сечения меняется по длине элемента. В слу чае линейного изменения площади величину А _в формулах
(12.8) — (12.10) |
можно заменить средней площадью A — (Ai-\-Aj)/2. |
Это выражение |
сразу же получается после замены А на А — |
= NiAi-\-NjAj и вычисления интеграла. Подобное выражение мо жет быть использовано для температуры, если она меняется ли нейно по длине. Нелинейные изменения учитываются с помощью интерполяционных полиномов, обсуждаемых в гл. 13. На следую щем примере показано, как видоизменить определяющие элемент соотношения, чтобы они соответствовали линейному изменению площади элемента.
Пример
108. Нужно вывести и решить систему линейных уравнений для узловых перемещений в конусообразной детали конструкции, один конец которой жестко закреплен, а другой подвержен действию на грузки в 42 000 Н. Площадь поперечного сечения меняется линейно от 12 см2 на левом конце до 6 см2 на правом. Кроме того, деталь конструкции испытывает тепловое расширение вследствие повы шения ее температуры на 20° равномерно по всей длине а = 7 Х XIО -6 1/°С. Для аппроксимации рассматриваемой части конструк ции следует использовать три элемента длиной 30 см каждый.
Площадь поперечного сечения в узловых точках имеет значе ния Ах= 12 см2, Л2= 10 см2, Л3= 8 с м 2 и Л4 = 6 см2. Первые три эле мента свободны и от объемных, и от поверхностных нагрузок, по этому матрицы этих элементов и векторы нагрузки определяются соответственно соотношениями (12.8) и (12.9).
К задаче 108.
Для первого элемента имеем
[£">[ = |
АЕ Г 1 |
— П |
_ 11-6,7* 1 0 е |
Г 1 |
— 1 |
|
L [ - 1 |
1J |
|
30 |
[ - 1 |
1 |
|
|
|
30 |
||||
|
= |
Г |
2,46 |
-2 ,4 6 ] |
|
|
|
|
[-2 ,4 6 |
2,46j |
|
{/(»} = а А Е (ЛТ) J— JJ = 7 • 10-«• 6,7 • 10е • 11 • 20 Jj J = j‘ •10318
10318
Уравнения, определяющие этот элемент, имеют вид
10еf 2'46 |
—2,46‘ |
£/П [— |
103181 |
[ - 2,46 |
2,46 |
t / j — 1 |
10318J* |
Первый и второй элементы различаются только размером пло щади поперечного сечения. Для второго элемента средняя пло щадь равна 9 см2. Определяющие уравнения для второго элемен та записываются как
Г 2,01 |
—2,01] р И = |
84421 |
[ —2,01 |
2,01 J [t/,1 |
\ 8442J |
Третьему элементу соответствуют уравнения, полученные с по мощью соотношений (12.8) и (12.9), а также (12.12), поскольку этот элемент нагружен на конце.
Г 1,56 — 1,56] f£/|| |
f—65661 + 7000-6 |
|
—65661 |
|||||||
[ - 1 , 5 6 |
1,56J \t/J |
[ |
6566J |
|
|
|
|
. 48566/ |
||
Соотношения включения: |
|
|
|
|
|
|
|
|||
для |
первого элемента: |
i — 1, |
/ = |
2, |
|
|||||
для |
второго |
элемента: |
i —2, |
/ = |
3, |
|
||||
для |
третьего |
элемента: t'= 3, |
у = 4 . |
|
||||||
Суммируя уравнения, определяющие элементы, получаем |
||||||||||
“ 2,46 |
—2,46 |
|
0 |
0 |
' |
\иА |
— 10318 |
|||
—2,46 |
|
4,47 |
—2,01 |
0 |
|
и, |
|
1876 |
||
0 |
—2,01 |
|
3,57 |
— 1,56 |
|
'— ' |
1876 |
|||
0 |
|
0 |
— 1,56 |
1,56 |
ы |
|
|
48566 |
Первый узел расположен в неподвижно закрепленной точке, поэтому Ui= 0, и приведенная система уравнений должна быть изменена с тем, чтобы учесть это граничное условие. В результате
имеем
“2,46 |
0 |
0 |
0 |
' [</.1 |
0 |
0 |
4,47 |
—2,01 |
0 |
V* |
1876 |
0 |
—2,01 |
3,57 |
— 1,56 |
|
' — ' 1876 |
0 |
0 |
— 1,56 |
1,56 |
|
48566 |
Приведем решение этой системы:
{£/}г=[0, 0,0207, 0,0450, 0,0753], см.
Теоретическое решение этой задачи получается путем интегриро вания деформации по длине. После выполнения этой процедуры получаем следующие значения для узловых перемещений:
t/1= 0 ,0 см, |
t/3= 0,046 |
см, |
U2—0,021 см, |
С/4 = 0,078 |
см. |
Перемещения, определенные методом конечных элементов, хоро шо согласуются с торетическими значениями. Еще более точные значения были бы получены при использовании элементов мень ших размеров.
12.1.1. Напряжения в элементах
Определение напряжений является важной частью решения большинства задач теории упругости, потому что эти величины используются инженерами для расчета различных элементов кон струкций. Результанты элемента, связанные с напряжениями, мо гут быть определены, как только вычислены деформации„внутри элемента. Для одномерной задачи деформация е** дается форму лой (12.6). Нормальное напряжение получается из закона Гука в форме (12.3).
Так как производные постоянны по элементу, деформация вну три отдельного элемента не меняется, что влечет в свою очередь в соответствии с законом Гука неизменность внутри элемента на пряжения. Узловые значения ахх могут быть рассчитаны с по мощью теории согласованных результантов элементов, представлен ной в гл. 6. Это делается аналогично тому, как было описано ра нее. Компоненты тензора напряжений являются результантами элемента. Теория согласованных результантов элементов может быть использована также для определения узловых значений ком понент тензора деформаций.
Пример
109. Для детали конструкции, рассмотренной в предыдущем примере, нужно рассчитать узловые значения ахх, используя тео рию согласованных результантов элементов.
Запишем вычисленные ранее узловые перемещения
|
{tf}r = [0 , |
0,0207, |
0,0450, |
0,0753]. |
|
||||
Определим теперь деформацию элементов: |
|
|
|||||||
первый элемент: |
ехх= - ^ - ( —U1 + U2) = |
° ’^ Q7*= 0 ,00069, |
|||||||
|
|
|
|
|
|
|
30 |
|
|
и |
элемент: |
ехх= |
—(/о |
Uо |
= |
—0,0207 + |
0,0450 |
л лплл « |
|
второй |
------^ |
3 |
----■’— 3~ |
’------ |
=0,00081, |
||||
третий |
элемент: |
ех |
-U3 + U t |
_ |
-0 ,0 4 5 0 + |
0,0753 |
_ Л Л Л 1 Л 1 |
||
30 |
|
---------------- |
|
30 |
|
0 ,0 0 1 0 1 . |
Напряжения в элементах даются формулой
охх= Е г хх—аЕ (АТ) = 6 ,7 • ЮвеХх—
— 7-10~®>6,7- 10е -20= 6,7 - 10вехх— 938.
Подставляя значения ех*, получаем
первый элемент: ахх= 3685 Н/см2,
второй элемент: а ^ = 4 4 8 0 Н/см2,
третий элемент: а** = 5820 Н/см2.
Уравнения теории согласованных результантов для элементов имеют вид
где |
Озсх— вычисленный результант для конкретного элемента, а |
Oi |
и оj — узловые значения ахх соответственно в узлах i и /. Запи |
шем эти уравнения для каждого элемента отдельно: первый элемент:
второй элемент:
третий элемент:
Объединим эти уравнения, используя метод прямой жесткости:
|
1 |
|
0 |
1 |
|
CN |
|
о |
|
1 |
1 2 |
1 0 |
||
6 |
0 |
1 2 |
1 |
|
|
||||
|
о |
О |
|
1SJC |
|
1 |
|
<*1 |
1842,5 |
<*2 |
4087 |
^3 |
5159 |
а 4. |
2914,5 |
Эта система имеет следующее решение:
{сг}г =[3558, 3935, 5222, 6132], Н/см2.
Теоретические значения напряжения оХх в узлах получаются делением величины приложенной нагрузки на площадь поперечно го сечения в соответствующей узловой точке. Три множества зна чений Охх приведены в следующей таблице:
|
Теоретическое |
Метод конечных элементов |
||
Номер |
|
|
||
узла |
значение |
согласованное |
напряжение, постоянное |
|
Охх= Р /А , Н/см2 |
||||
|
напряжение, Н/см2 |
по элементу, Н/см2 |
||
1 |
3500 |
3558 |
|
|
2 |
4200 |
3935 |
Первый элемент 3685 |
|
3 |
5250 |
5222 |
Второй элемент 4489 |
|
4 |
7000 |
6132 |
Третий элемент 5829 |
Значения ахх» вычисленные по теории согласованных результан тов, определенно лучше значений напряжения, постоянных по эле менту, но они все же еще недостаточно близки к теоретическим значениям. Дальнейшее улучшение значений ахх может быть до стигнуто путем применения элементов меньших размеров.
12.2. Двумерные задачи теории упругости
Двумерные задачи теории упругости намного сложнее одно мерных, поскольку в случаях плоского напряженного или плоско го деформированного состояния может иметь место анизотропия материала. Каждому из этих двух состояний соответствует своя
матрица упругих характеристик [D].
В плоских задачах теории упругости применим треугольный симплекс-элемент с шестью компонентами узловых перемещений
(фиг. 12.1). Перемещения |
и н и |
внутри |
элемента даются зависи |
|||||
мостью |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
UU-i |
|
|
М |
Г Nt |
0 |
| N, |
О |
Nk |
U2j- 1 |
(12.14) |
|
Ы |
[О |
А7, |
О |
N} |
О |
U2j |
||
|
U*k-i
V2k
Обозначения узловых перемещений показаны на фиг. 12.1. Функ ции формы, входящие в соотношения (12.14), определены в (3.10).
Будем считать, что рассматриваемая область располагается в плоскости ху, и введем следующие компоненты напряжения и де
формации! {о}‘Г== [Охх> Оуу> Тосу] И { б } Т = [бхх> 6уу> У ху] • Д Л Я ПЛО-
Фиг. 12.1. Компоненты перемещения для двумерного симплекс-элемента.
ского напряженного состояния, встречающегося во многих тонких
телах, имеем |
azz= T zx= T zl/= 0 . Компоненты тензора деформации |
уyZ и уХг тоже |
равны нулю, но ezz отлична от нуля и может быть |
получена из закона Гука, после того как определены {а} и {е}. Говорят, что плоское деформированное состояние имеет место, ког да компоненты деформации в направлении оси z равны нулю (ezz=Yxz=yyz=0). Компоненты тензора напряжений т2у и т2Х так же равны нулю при плоской деформации, но а22 отлична от нуля и вычисляется с помощью закона Гука после того, как определены {а} и {е}.
Соотношения связи между деформациями и перемещениями в двумерном случае имеют вид
|
ди |
|
|
да |
^ХУ~~ |
ди |
|
да |
|
|
|
дх |
|
г уу |
"ду ’ |
ду |
|
дх |
|
||
или с учетом (12.14) |
|
|
|
|
|
|
|
|
|
|
Чс* |
|
bi 0 |
bj |
0 |
bk |
0 ' |
v* |
|
|
|
|
1 |
0 |
|
0 |
|
0 |
|
Uv- 1 |
|
|
г уу |
2А |
|
Ч |
|
ч |
|
4 |
и* |
(12.15) |
|
Еху, |
|
|
Ч |
CJ |
bj |
4 |
bk_ |
i |
|
|
|
|
|
|
|
|
|
|
f V |
|
|
|
|
|
|
|
|
|
|
а.2к |
|
|
Соотношения |
(12.15) |
|
определяют |
матрицу |
градиентов |
[В], |
||||
так как {е} = [£ ]{[/} . Теперь есть |
почти |
все необходимое для |
вы |
вода уравнений, определяющих элемент. Осталось только записать матрицу упругих характеристик [Z?] и вектор начальной деформа ции {е0}. В случае плоского напряженного состояния имеем
~1 р О
|
[D] = |
|
р |
1 |
0 |
, |
|
|
|
О |
0 |
(1— р)/2 |
|
|
|
(е0)= а Д 7 ’ |
Т |
|
|
|
|
|
1 |
|
|
||
|
|
|
|
О |
|
|
В случае плоской деформации |
|
|
|
|
||
\D] = |
Е( 1 -р ) |
'1 |
р /(1 - р ) |
О |
||
р/(1— р) |
1 |
|
|
О |
||
(1+ р )(1 -2р ) |
|
|
||||
|
|
О |
О |
|
|
(1— 2р)/2 (1— р) |
и
Т
(е0} = ( 1 + р )а Д Г 1
О
(12.16)
(12.17)
(12.18)
(12.19)
Формулы (12.16) — (12.19) соответствуют изотропному материалу с модулем упругости Е и коэффициентом Пуассона р.
Интегралы, на основе которых составляются уравнения, опреде ляющие элемент, легко вычисляются, поскольку матрицы [В] и
[Щ содержат только константы. Вычислим объемный интеграл, представляющий матрицу жесткости:
[£M]=j* [BY [D] [В] d V = [B f ID] [В] j* dV,
V V
[№]={B\T[D\[B\tA. (12.20)
Здесь t — толщина элемента, А — его площадь. Общее выражение для матричного произведения [B]T[D][B] не приведено из-за его громоздкой записи. Обычно поступают так: определяют числовые значения коэффициентов [В] и [D], а затем ЭВМ выполняет ука занное перемножение матриц. Интеграл, связанный с тепловым расширением, имеет вид
[B]T[D] {е0} dV = - [ B ] T[D]E0\ tA. |
(12.21) |
и
Матричное произведение в формуле (12.21) нетрудно составить. Для случая плоских напряжений получаем
|
|
ЬЛ |
|
|
Ci |
[B\T[D][e0}tA |
аЕЦЬТ). Ьг |
|
2 0 |
( 12.22) |
|
|
Ц |
bk
Объемный интеграл от объемных сил аналогичен интегралу
j[N]T{Q}dV, который был получен при рассмотрении задач тео
рии поля. Основное отличие заключается в том, что теперь матрица ;[W]r состоит из двух столбцов, так как имеются две объемные си лы. Подставляя [УУ]Г и применяя /^координаты, получаем
■w« |
0 |
" |
|
Ц |
|
01 |
Nt |
|
|
У |
|
0 |
f |
dy— L |
зе |
(12.23) |
|
0 |
Nj |
Ы |
3 |
У |
|
Nk |
0 |
|
|
эе |
|
0 |
Nk |
|
У |
|
Интеграл от поверхностных нагрузок также аналогичен поверх ностному интегралу в задачах теории поля. Рассматривая отдель но каждую из сторон элемента, можно записать три различных значения этого интеграла. Предположив, что на стороне между узлами i и / действуют равномерно распределенные нагрузки ин