книги / Применение метода конечных элементов к расчету конструкций

ГЕОМЕТРИЧЕСКАЯ МАТРИЦА ЖЕСТКОСТИ ОБОЛОЧКИ. Рассмотрим построение геометрической матрицы жесткости на примере изотропной упругой оболочки произвольной формы, представ­ ляемой ансамблем плоских треугольных КЭ постоянной толщи­ ны (рис. 6.5.2).

Рис. 6.5.2. Конечноэлементная модель тонкой оболочки*

Запишем зависимость вектора перемещений и(х^) в любой точке КЭ от вектора обобщенного перемещения узлов v треу­ гольного КЭ:

и(х^ )=G(x-1)v.

(6.5.7)

Рис. 6.5.3. Плоский изгибаемый треугольный КЭ.

Выделим отдельно мембранные и изгибные перемещения. Обоз­ начим через Up и и^, соответственно, векторы обобщенвдх

мембранных и изгибных перемещений. Тогда uT=[Up u^].

Аналогичное обозначение введем и для матрицы функций фор­

мы: GT=[Gp Gb ],

Компонентами вектора перемещения v являются перемеще­ ния в плоскости КЭ и перемещения изгиба (рис. 6.5.3):

v T- [Ui

w.

l )i

X

Вектор деформаций КЭ равен:

V			------------------1	р	м
ех2				V,2
е= V	x 2			" , 2 +\ l
V				W, l l
х				W,22
У Г
*	1. 2	J	•	2W,12
V	V

(<р 2 )i		(»> 3 )±
X		X
1	(w	, ) 2/2
_________
	f	" L
	(wf 2 ) 2/ 2
+		О
		Г
		0
		0

-0

...].

1

' С	' 4 '		(6.5.8)
	- Л . 1.
г,		поверхности;		0 _
Здесь ц, v, w - перемещения срединной				£р
вектор мембранной деформации (первые 3 компоненты);
вектор деформации растяжения за счет		изгиба;	-	вектор
изгибной деформации.
Потенциальная энергия деформации КЭ равна:
n=T	[}jy<jdv,		( 6 . 5 . 9 )
где
	ег=Се,		(6.3»10^

Т	<r2	T		1	2	m		l	m		2	m	1 2	(6.5.11)
(Г— crl							X			X
X	X		X X										X X
Для изотропного		материала:
	г 1			Г		0				0			0	0
	Г			1		0				0			0	0
Е	0			0	(1-г )/2					о			0	0
с=	0			0		0			h3/12 rh3/l2					0
i-„2
	0			0		0			rh3/12				h3/l2	0
	- 0			0		0				0			0 (1-■r)h3/l
Е - модуль упругости,				г	- коэффициент Пуассона, h - тол­
щина плоского	изгибаемого					КЭ.			в выражение					(6.5.9). После
Подставим	(6.5.8),			(6.5.10)

интегрирования по толщине пластины и отбрасывания слагае­

мых, дающих вклад в энергию деформаций больших								прогибов
как малых более высокого порядка					(считается,		что	переме­
щения	малы), потенциальную			энергию П можно		представить в
виде:
П=	Eh				1-г
	.<U ,l*2+(V ,2>2+2l’U/lV,2+ - 7 - < V, r u ,2>2]dA+
2 ( 1 - 0 J
	Eh'				2 2 ^ +2vW 11W
	+ -------------- о-	I I M W	1 1 > 2 + (W				2 2 +
	24(l-r ) JJL		'11		'X1		•jLZ
	+2(l-r) (w^12)2JdA+-—			j|^<
				a lh(w/l^2+0' 2h<W,2)2+
				X	X
				A
	+T	. „hw	,w ^IdA,					(6.5.12)
		x V	,1	,^J
где A	- площадь КЭ; a -h,			a 9h,	т - 0h - интенсивности
		xA		x	x’Sc	усилия), т.е.
усилий в срединной плоскости (мембранные
силы на единицу длины.

С учетом зависимости (6.5.7) потенциальная энергия равна:

1 rt !I=1	г
	н>
2	L о

О

1

+

KbJ

1

ОО

о

£

ч

(6.5.13)

J

Здесь Кр - матрица мембранной жесткости КЭ:

R	----- т ГГГо		-G*	,+G	^ +P 9[G		-G*	9+G	,G * . ] +
P		JJ L Pf^ P/l		P r2	pf2	( Pfl	Pr2	P/2	P/1J
	l-l>
		(ep.24V l ) ( ep.2+ep , l ) ' h '							(б-5 Л 4 >
- матрица		изгибной жесткости КЭ:
	Eh		Г						ГП
КЬ	24(1-v2.) JJ [GbrllGb/ll+Gbr22Gb/22^{Gb fllGb,22+
	+Gb,22Gb,ll)+2<l-'>Gb,12Gb,12]dfl!								(6.5.15)
KJJ - матрица		геометрической жесткости КЭ:
			[axlhGb ,lGb, l+ffx2hGb ,2Gb ,2+
	+т	4	Jд						(6.5.16)
		X X	[Gb ^ Gt U +Gb,2Gi!U)]di-

УРАВНЕНИЕ УСТОЙЧИВОСТИ СТРУКТУРЫ. Потенциальная энер­ гия деформации структуры равна сумме потенциальных энер­ гий деформаций отдельных КЭ. Введем следующие обозначения для i-ro КЭ:

о

К. = 1

О

«У

Потенциальная энергия деформаций структуры равна:

П= — v (К+К )v,

(6.5.17)

2 и

где	m
m

к = l Li Ki Lr	ки - I LI KMiLi :
i=l	i=l
- матрица направляющих	косинусов для i-ro КЭ; m -
число КЭ в структуре.

Работа, совершаемая нагрузками f, приложенными в узлах

	структуры имеет вид:	1	т
	W			(6.5.18)
		= —	vf.
	®	2

При простом однопараметрическом .нагружении узловые нагрузки одновременно изменяются пропорционально одному параметру Л. Это означает, что для произвольного уровня нагружения имеет место формула:

f=XfQ,

(6.5.19)

где fg - вектор начальных узловых нагрузок.

Матрица геометрической жесткости структуры

зависит

от уровня достигнутых напряжений или усилий при увеличе­ нии параметра Л, определяемых изменением конфигурации структуры. Учитывая, что эта зависимость является линей­ ной, запишем уравнения равновесия на основе выражения (6.5.17) и (6.5.18) в следующем матричном виде:

Kv+AKN (tf0)v=\fQ,

(6.5.20)

где Фо=Фо[*о) “ вектор напряжений

изменение параметра Л приводит к изменению изгибной жесткости структуры, которое обусловлено вторым слагаемым в левой части равенства (6.5.20). В случае оболочки это изменение обусловлено действием тангенциальных напряжений в ее срединной поверхности, в случае изгибаемого стержня - действием осевой растягивающей или сжимающей силы. При действии растягивающей силы изгибная жесткость стержня увеличивается, а при действии сжимающей силы изгибная жесткость уменьшается. Считаем, что эти напряжения или сила являются консервативными, т.е. совершаемая ими рабо­ та на любых кинематически возможных перемещениях зависит только от начальной и конечной конфигурации структуры.

Потере устойчивости структуры соответствует такое значе­ ние параметра А, при котором решение матричного уравнения (6.5.20) является не единственным. В этом случае справед­ ливо уравнение:

Kv+XKN (^0)v=0/	(6.5.21)
которое является уравнением устойчивости структуры.
Это уравнение аналогично уравнению (6.3.6)	для, нахож­

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

МЕТОД ОПРЕДЕЛЕНИЯ КРИТИЧЕСКОГО ПАРАМЕТРА НАГРУЗКИ X С ИСПОЛЬЗОВАНИЕМ ПРИЗНАКА ДАЛАМБЕРА СХОДИМОСТИ ЧИСЛОВОГО

РЯДА. В уравнении (6.5.20) вектор if/Q		является функцией
вектора сил f^. Рассмотрим линейное уравнение:
решение которого имеет	вид:
v 1=K-1f().		(6.5.22)
Вектор перемещений	используется	для вычисления	и
затем матрицы геометрической жесткости

Умножим обратную матрицу К ^ слева на обе части

уравнения (6.5.20). С учетом равенства (6.5.22) получим:

(I +A K " 1K N )V =A V 1 ,

откуда

V =(I+AK“1KN )“1AV 1.

Здесь I -единичная матрица.

Разложим (I+AK ^KJJ) * в матричный степенной ряд:

СО

(I+AK" 1KN )“1=1+ 5](“1)П”1АП"1(К“1KN )П"1.

п=2

С учетом этого ряда для вектора перемещений v имеем:

00

v=Av1+[^ ^(-l)n"1An (K"1KN )n"1jv1.

(6.5.23)

n=2

Введем следующие обозначения: vn=-K *K N V II-1 (n=2,3,...) .

Векторы v2 , v 3 f ..,,vn являются приращениями узловых

перемещений структуры, обусловленные матрицей геометриче­ ской жесткости KN » Они зависят только от начальной наг­

рузки £ и	не	зависят от уровня	нагружения,	определяемого
параметром	нагрузки А .		ряд (6.5.23)	запишется в
С учетом		этих обозначений
виде:		00
		V = £*nvn .		(6.5.24)
		п=1

При приближении параметра нагрузки А к критическому значению Акр в процессе увеличения уровня нагружения ряд

(6.5.24) будет расходиться, т.е. вектор узловых перемеще­ ний v будет неогранниченно возрастать, что соответствует состоянию потери устойчивости структуры в смысле Эйлера. Абсолютное значение критической нагрузки тогда равно:

fKp“XKPf0*

Для численного определения Лкр воспользуемся признаком Даламбера сходимости числового ряда. Обозначим:

.п+1„

АиV.п+1

		An llv_
		п
где Ivll -норма вектора узловых перемещений структуры.
При	7)<1 ряд (6.5.24)	будет сходящимся. При 7]=1 ряд
будет	расходиться. Условие расходимости ряда позволяет
вычислить:
	A, “llv ll/llv
	ът	П /	п+1
	кр

Форма потери устойчивости находится в результате вы­ числения суммы ряда:

	п+1
v -	= £A _._V	k*
кр	кр
	k=l

Описанный метод применим и к задаче нелинейного деформи­ рования структуры.

Г л а в а 7

РЕАЛИЗАЦИЯ МКЭ НА ЭВМ

7.1. ПРОГРАММНЫЕ СИСТЕМЫ

Развитие метода конечных элементов обусловлено взаимо­ связью трех факторов: наличием высокопроизводительной вы­ числительной техники; разработкой математических моделей исследуемых явлений, адекватных реальным процессам с достаточной степенью точности; особенностями самого мето­ да [76].

Первые программные комплексы, в которых реализован ме­ тод конечных элементов, были разработаны в б0-х годах [200]. К ним относятся STRUDL-II, SAP-IY, NONSAP, ASKA, NASTRAN, SESAM-69 и другие. Появлению этих универсальных программных систем .я силу особенностей метода конечных элементов предшествовало создание высокопроизводительных электронно-вычислительных машин, таких, например, как IBM-370. Начиная с конца 70-х годов в СССР появилось не­ сколько десятков программных комплексов для разных ЭВМ, в которых был реализован МКЭ. К их числу относятся МИРАЖ [39], МОРЕ [78], ПРАСАК[62], КАСКАД-2 [86], ПРОЧНОСТЬ-75 [55], СИСТЕМА-4 [64],СИПРАМАК [30], МКЭ/120 [27], МАРС [32], АПЖБК [82], ПАРДОКС [23], ПАРСЕК [23], ЛИРА [41], СУМРАК [17], СПРИНТ [115], ПРОКРУСТ [16], ПАРУС [28],ПОЛИФЭМ [46] и ряд других программ.

В США и ряде других стран дальнейшее развитие МКЭ и необходимость в проведении расчетов конструкций на проч­ ность также способствовали дальнейшему развитию уже соз­ данных программных комплексов и разработке новых. .Были разработаны сотни программных комплексов, предназначенных для приближенного решения самых разнообразных задач не

только из области механики деформируемого твердого	тела,
но и из таких областей как гидродинамика, акустика,	элек­

тротехника и т.д. Наибольшее распространение из них полу­

чили [193]: ABAQUS,	ADINA, ASKA/DYNAN, ANSYS, FINEL,
LARSTRAN, MSC/NASTRAN,	EUFEMI, ROSALIE, HERCULE, MODULEF,
SAP-7.
Отметим, что разработка программных комплексов являет­
ся дорогостоящим делом. Поэтому, как правило,		организации
и фирмы - собственники разработанных программ,		рассматри­

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

диные для их деятельности программы расчета.

У каждой программы есть свои сильные и слабые стороны при расчете конкретной конструкции. Выбор программы рас­ чета зависит от подготовленности пользователя в своей научной области, типа решаемой задачи, типа доступной ЭВМ, размерности задачи и других факторов. К критериям, помогающим сделать выбор, следует отнести следующие фак­ торы:

- программа широко используется; - в программе использованы новейшие научные достижения;

-программа коммерчески вполне доступна;

-имеется подробная и понятная документация.

.Ознакомление с программной документацией и доступной литературой с описанием программы и ее элементов позво­ ляют сделать окончательный вывод о целесообразности выбо­ ра программного комплекса.

Для МКЭ характерны особенности, которые следует учиты­ вать при выборе или разработке программы расчета. Такими особенностями являются большие объемы исходных данных, промежуточных и окончательных результатов расчета. Поэто­ му расчет по МКЭ состоит из трех основных этапов:

-разработка расчетной конечноэлементной схемы и подго­ товка исходных данных;

-проведение самого расчета;

-обработка результатов расчета.

На рис. 7.1.1. приводится одна из возможных схем орга­ низации расчета по МКЭ. Каждый этап является самостоя­ тельной задачей. На первом этапе самое существенное за­ ключается в создании начальной конечноэлементной расчет­ ной модели, исходя из инженерной интуиции о поведении конструкции. Впоследствии эта модель может корректиро­ ваться на основе анализа результатов расчета. Корректи­ ровка модели может выполняться и программным путем, если такая возможность реализована в используемом программном комплексе. Подготовка исходных данных осуществляется, как правило, с помощью программ-генераторов сеток конечных элементов, образующих блок подготовки данных. Подробное описание возможностей таких программ блока обработки ре­ зультатов расчета (этап 3) будет приведено в разделе 7.3.

Проведение расчета (этап 2) осуществляется расчетным блоком, в котором используется тот или иной алгоритм расчета методом конечных элементов. Как правило, расчет­ ный блок состоит из ряда программных модулей, каждый из которых выполняется на определенном шаге алгоритма, в простейшем случае программной реализации МКЭ для линейной

статической краевой	задачи теории	упругости расчетный
блок содержит следующую последовательность шагов:
- ввод исходных данных (например,		подготовленных прог­
раммой-генератором в	отдельном файле);

автоматическая подготовка исходных данных

диагностика качества сетки КЭ с использованием графических средств

необходимо качество сетки улучшить

корректировка сетки КЭ

диагностика необходимости уменьшения ширины ленты матрицы жесткости

необходимо

перенумерация узлов сетки КЭ

расчетный блок

анализ результатов и диагностика точности

точность удовлет.

вывод результатов i) табличном и (или) графиче<жом виде

конец вычислений

Рис. 7.1.1. Схема организации расчета по МКЭ.

-вычисление матриц жесткостей КЭ;

-формирование матрицы жесткости структуры;

-формирование вектора нагрузок;

-решение системы линейных алгебраических уравнений;

-вычисление перемещений узлов ансамбля конечных элемен­ тов, деформаций и напряжений в конечных элементах.

Пример организации расчетов показан на рис. 7.1.2.