Спецкурс численные методы.

											81
Для первой функции:
a1 *0 b1	1,
a1 *l b1	0,
откуда b 1, a				1			и , следовательно,		(x) 1	l	x 1.
1		1				l			1
Для второй функции
a2 *0 b2	0,
a2 *l b2	1,
откуда b 0, a		2	1		l	и , следовательно,		2	(x) x l .
2

Пример 2. Дано: одномерная область, порядок дифференциального уравнения s 4. В этом случае число степеней свободы КЭ N s 4 . Конечный элемент - отрезок прямой длиной l. Задавая на каждом конце КЭ по два параметра - значению искомой функции и её производной на границе элементов, мы полностью определим приближённое решение на носителе. Выберем в качестве координатных функций многочлены степени p s-1 4-1 3, то-есть, кубические многочлены.

Общий вид кубического многочлена (x) a	i	x3	b	i	x2	c	i	x	d	, его производная -
i									i

i (x) 3ai x2 2bi x ci x.

Первый из них удовлетворяет условиям

1(0) 1, 1 (0) 0, 1(l) 0, 1 (0) 0 ,

второй –

2(0) 0, 2 (0) 1, 2(0) 0, 2 (0) 0,

третий –

3(0) 0, 3 (0) 0, 3(l) 1, 3 (0) 0,

четвёртый –

4(0) 0, 4 (0) 0, 4(l) 0, 4 (l) 1 .

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

a1*03 b1*02 c1*0 d1 1, 3a1*02 2b1*0 c1 0, a1*l3 b1*l2 c1*l d1 0, 3a1*l2 2b1*l c1 0 ;

a2*03 b2*02 c2*0 d2 0, 3a2*02 2b2*0 c2 1,

a2*l3 b2*l2 c2*l d2 0,

82

3a2*l2 2b2*l c2 0;

a3*03 b3*02 c3*0 d3 0, 3a3*02 2b3*0 c3 0,

a3*l3 b3*l2 c3*l d3 1, 3a3*l2 2b3*l c3 0;

a4*03 b4*02 c4*0 d4 0, 3a4*02 2b4*0 c4 0,

a4*l3 b4*l2 c4*l d4 0, 3a4*l2 2b4*l c4 1.

Вектора решений этих систем представим в виде таблицы (табл. 7.1).

				Таблица 7.1.
Номер	Коэффициенты многочлена – i-ой координатной функции
координатной	ai	bi	ci	di
функции
1	2/l3	-3/l2	0	1
2	1/l2	-2/l	1	0
3	-2/l3	3/l2	0	0
4	1/l2	-1/l	0	0

Введём удобное для анализа свойств координатных функций относительное независимое переменное xl . Тогда четыре координатные функции этого КЭ примут вид:

1( ) 2 3 3 2 1;2( ) l (1 )2 ;3( ) 2(3 2 );4( ) l 2(1 ).

Пример 3. Дано: двумерная область, порядок дифференциального уравнения по

каждой координате	sx sy s 2. Тогда число степеней свободы	КЭ N 2s 2*2 4 .
Конечный элемент	– прямоугольник x 0,l , y 0,h . Задавая на	каждом узле КЭ по

одному параметру - значению искомой функции, мы полностью определим приближённое решение на носителе. Выберем в качестве координатных функций многочлены степени p s-1 2-1 1 по каждой координате (то-есть, функции, линейные по каждой координате). Общий вид таких функций:

i (x, y) aixy bxi ci y di .

83

Первая из них удовлетворяет условиям

1(0,0) 1; 1(0,1) 0; 1(1,1) 0; 1(1,0) 0;

вторая -

2(0,0) 0; 2(0,1) 1; 2(1,1) 0; 2(1,0) 0;

третья –

3(0,0) 0; 3(0,1) 0; 3(1,1) 1; 3(1,0) 0;

четвёртая –

4(0,0) 1; 4(0,1) 0; 4(1,1) 0; 4(1,0) 1.

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

a1 *0*0 b1 *0 c1 *0 d1 1, a1 *0*h b1 *0 c1 *h d1 0, a1 *l*h b1 *l c1 *h d1 0, a1 *l*0 b1 *l c1 *0 d1 0;

a2 *0*0 b2 *0 c2 *0 d2 0, a2 *0*h b2 *0 c2 *h d2 1, a2 *l*h b2 *l c2 *h d2 0, a2 *l*0 b2 *l c2 *0 d2 0;

a3 *0*0 b3 *0 c3 *0 d3 0, a3 *0*h b3 *0 c3 *h d3 0, a3 *l*h b3 *l c3 *h d3 1, a3 *l*0 b3 *l c3 *0 d3 0;

a4 *0*0 b4 *0 c4 *0 d4 0, a4 *0*h b4 *0 c4 *h d4 0, a4 *l*h b4 *l c4 *h d4 0, a4 *l*0 b4 *l c4 *0 d4 1.

Вектора решений этих систем представим в виде таблицы (табл. 7.1).

Таблица 7.1.

84

Номер	Коэффициенты многочлена – i-ой координатной функции
координатной	ai	bi	ci	di
функции
1	1/(lh)	-1/l	-1/h	1
2	-1/(lh)	0	1/h	0
3	1/(lh)	0	0	0
4	-1/(lh)	1/l	0	0

Введём удобные для анализа свойств координатных функций относительные независимые переменные xl , yh . Тогда четыре координатные функции этого КЭ примут вид:

1( , ) 1;

2( , ) (1 ) ;

3( , ) ;

4( , ) (1 ) .

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

Существенно, что в наших задачах количество главных и естественных краевых условий одинаково и это энергетически сопряжённые величины – интеграл от произведения естественного условия на вариацию главного условия равен вариации работы обобщённого усилия (естественного условия) на соответствующем возможном обобщённом перемещении (главном условии). В процессе расчёта подбираются такие значения главных краевых условий в узлах (одинаковые для всех узлов КЭ, которые совмещены с данным узлом системы), чтобы в этих же узлах выполнялись (в приближённом смысле) естественные краевые условия. В задачах строительной механики это означает, что узлы должны находиться в равновесии, при этом условия равновесия могут записываться в трёх эквивалентных вариантах:

либо в виде равенства нулю обобщённой силы, действующей на узел, либо в виде равенства нулю работы обобщённой силы на соответствующем

обобщённом перемещении, либо как условие минимума функционала энергии.

Если среди краевых условий есть главные и естественные, то, как мы выяснили, главные условия должны выполняться точно, а естественные – приближённо. В прямых методах это означает, что главные условия должны выполняться уже во время конструирования координатной системы, а приближённое соблюдение естественных уловий должно обеспечиваться в результате расчёта. Исходя из этого, условия сопряжения подобластей в узлах обеспечивается:

-строгим равенством величин, задающих главные условия (в строительной механике – обобщённых перемещений нормали и её поворотов для балок и плит и перемещений для массивов);

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

85

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

Строгое равенство главных условий в налагаемых друг на друга узлах смежных элементов (можно сказать – в совмещаемых узлах смежных элементов) обеспечивается в алгоритме (и программе) тем, что неизвестными являются приближённые значения

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

Рассмотрим, какой вид приобретают коэффициенты системы уравнений методов Ритца или Галёркина, когда эти методы сочетаются с МКЭ. Согласно традиции, установившейся в литературе по МКЭ, будем обозначать эти коэффициенты как kij , если

речь идёт о всей системе, и kij(e) , если речь идёт о КЭ с номером e (вместо rij по ЧМ-6). Вначале выберем КЭ, занимающий всю рассматриваемую область.

Введём обозначения: e - номер выбранного КЭ; Ie - множество номеров узловых параметров этого КЭ (нумерация в пределах КЭ); Be - носитель всех координатных функций, связанных с выбранным КЭ; B - вся область (вначале Be B ). Согласно (6.30), обобщая эту

формулу, в случае необходимости, на случай пространств большей размерности, и обозначая аргумент аналогично тому, как это сделано в пояснениях к формуле (7.1), получим

kij(e) L j (

x

) i

(

x

) d(

x

) L j (

x

) i(

x

) d(

x

).

(7.2)

B

Be

Этот интеграл

в случае неодномерной области – кратный,

d(

x

) dx1dx2 или

d(x) dx1dx2dx3 . При вычислении этого интеграла может быть выполнено понижение

порядка дифференциального оператора вдвое путём применения интегрирования по частям или его обобщения для неодномерной области (формул Грина, Стокса или ГауссаОстроградского) или с учётом физического смысла соотношений метода Ритца..

Количество элементов множества Ie		- это количество узловых параметров КЭ, число
его степеней свободы Ne	. Наборы ij Ne	при каждом фиксированном значении j const
	i 1

Ne
можно считать j - ми столбцами, а наборы ij j 1	при каждом фиксированном значении

i const i- ми строками некоторой матрицы K(e) kij(e) . Эта матрица называется матрицей

жёсткости конечного элемента № e - МЖ КЭ № e. Если оператор задачи симметричен, то эта матрица также симметрична. Если КЭ представляет собой неклассическую область, вычисление интегралов (7.2) может представлять собой нетривиальную и уж во всяком случае трудоёмкую задачу. Все промышленные программы МКЭ имеют библиотеки конечных элементов с формулами для определения элементов МЖ КЭ.

В ряде случаев различные КЭ одной конечно-элементной сетки могут быть сориентированы по-разному. В этих случаях целесообразно вводить на каждом КЭ свою систему координат. Такая система координат называется местной (локальной) системой координат данного КЭ.

Если ввести общую (глобальную) систему координат, то взаимное положение этой системы и каждой из местных систем может быть охарактеризовано углом взаимного

86

поворота e . Тогда компоненты векторов Z(e) , Q(e) , заданные в местной системе, могут быть преобразованы к общей системе по формуле, повторяющей преобразование координат

Ze(g) T(e)Ze(l) ,			(7.4)
Q(g) T(e)Ql		,	(7.5)
e	e

где T(e) - матрица преобразования координат (матрица поворота) от местной системы к глобальной; Ze(l) , Ze(g) - узловые параметры КЭ, соответствующие главным граничным

условиям; Q(l),	Q(g)	- узловые параметры КЭ,	соответствующие естественным граничным
e	e
условиям,
В аналитической геометрии доказывается, что T(e)				- ортогональная матрица, то-есть,
она обратима и её обратная матрица равна её транспонированной матрице
[T(e) ] 1 [T(e)]T .					(7.6)
Тогда из (7.4) с учётом (7.6) можно записать
Ze(l) [T(e)] 1Ze(g) [T(e) ]T Ze(g) .					(7.7)
Подставим (7.7) в (7.3), после чего - преобразованное выражение (7.3) в (7.5)
Q(g) T(e)Q(l)		T(e)K(e)Z(l) T(e)K(e)[T(e)]T	Z(g) K(g)Z(g) ,		(7.8)
e	e	e	e	e
где K(g)	T(e)K(e)[T(e) ]T - МЖ КЭ в глобальной системе координат.

В МКЭ принято определять естественные граничные условия КЭ, приведенные к его узлам (узловые обобщённые силы КЭ в случае задачи строительной механики)

P (( f , ),( f ,			),...,( f ,	))T .	(7.9)
e	1	2		Ne

Теперь перейдём к случаю сетки КЭ с произвольным количеством элементов. Узлы, принадлежащие какому-либо КЭ либо совпадающие с его узлами, в МКЭ называют инцидентными ему. Соответствие между номерами КЭ и номерами инцидентных им узлов системы задаётся расчётчиком в исходных данных в виде специальной таблицы инцидентций.

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

		ne	qn
			(l)			...,	(7.10)
y(x) ... ze,nq (x)e,nq
		e Em n 1	q 1

где e- номер КЭ; Em - множество номеров КЭ, которым инцидентен узел m системы; n-

номер узла КЭ в его местной нумерации; ne - количество узлов в e-ом КЭ; q- номер

87

узлового параметра данного узла; qn -количество узловых параметров в n-ом узле; ze(l,nq) - q-

ый узловой параметр n- го узла; (x)e,nq - координатная функция, соответствующая этому

параметру; многоточием показаны нулевые слагаемые (связанные с нулевыми узловыми параметрами).

Пусть теперь в некотором узле системы с глобальным номером m задан своими координатами в глобальной системе координат вектор узловых параметров Zm(g) . По аналогии

с (7.10), вспоминая понятие скалярного произведения векторов из линейной алгебры, запишем (с учётом (7.7)) соответствующий вектор аппроксимаций

	Y(	x	) Te,n(m) (	x	)Ze(,ln)(m) Te,n(m)(	x	)T(e)T Zm(g)	e Em ,	(7.11)
где	Te,n(m) - транспонированный вектор							координатных функций,	соответствующих

заданному вектору узловых параметров; n(m)- номер узла КЭ, накладываемого на m-ый узел системы.

Учитывая, что за пределами носителя координатные функции тождественно равны нулю, (7.11) можно записать в виде

Y(	x	) Te,n(m) (	x	)T(e)T Zm(g) .	(7.12)
		e Em

По аналогии с формулами (6.30), (6.38) и с учётом (7.8) здесь следует иметь в виду, запишем блоки матрицы коэффициентов МКЭ в форме методов Ритца или БубноваГалёркина (опять согласно традиции МКЭ заменим обозначение коэффициентов на k ; блок

K

pq

k

1 qn 1, 1 qn 1

, int(

),

int(

))

1, 1

1

qn

1

qn

T

(l)

K

)

)

(

Z

))dx

(7.13)

m

(LY ,

((LY

e,

e,

m

m

Be Em

( ((L e,n(m) )T(e)T E)T ( e, T(e)T E))dx

Be Em

T(e)( (L e,n(m))T e, dx)T(e)T Ke(,gm)

e Em	B	e Em
	e

Здесь E	m	E	m	E ;	K(g)	- отдельные блоки МЖ КЭ e, расположенные на
					e,m

пересечении её m-ой блочной строки и -го блочного столбца; Nm , где

Nm - множество

номеров узлов системы, инцидентных тем же КЭ, что и узел m, Nm .

Для столбцов свободных членов

аналогично (7.9) получаем выражение для m-го

блочного подвектора

Pm f (

x

) m (

x

)dx

f (

x

) e,m (

x

)dx

Pe,m .

(7.14)

E

e Em E

e Em

m

m

88

Для матрицы коэффициентов всей системы уравнений для определения узловых параметров рассчитываемого объекта и её столбца свободных членов естественно получаем,

используя (7.13) и (7.14)

K	Km	,	(7.15)
P Pm .			(7.16)
Упомянутая система уравнений в матричной записи имеет вид
KZ P.			(7.17)

Это разрешающая система уравнений МКЭ.

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

Решая систему (7.17), получаем приближённые значения искомой функции в узлах системы

	,	(7.18)
y(xm ) zi

m 1

где i qn 1, m- номер узла.

n 1

Её производные порядка до qn 1 включительно (здесь qn - количество узловых параметров n- го узла системы) определяются по формуле (7.18) при

m 1
i qn p 1,	(7.19)
n 1

где p qn 1 - порядок разыскиваемой производной.

Обычно необходимы не только значения искомой функции в отдельных узлах системы, получаемые в результате решения системы (7.17), но также и значения функции и её производных (или их линейных комбинаций) между узлами, то-есть, внутри КЭ. В этих случаях для каждого КЭ строятся аппроксимации искомой вектор-функции (7.11).

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

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

89

Следует подчеркнуть, что предварительная (априорная) оценка погрешности практически невозможна, и единственная возможность получить приближённое решение с погрешностью, не превосходящей заданную – организовать последовательность процессов решения задачи при сгущении узлов и измельчении КЭ от расчёта к расчёту. При грамотной подготовке исходных данных последовательность полученных приближённых решений должна сходится к точному решению, а разность двух решений может послужить для оценки погрешности (апостериорной), если, например, сгущение узлов осуществлять только путём добавления новых узлов и сравнивать вектора значений неизвестной функции в тех узлах, которые есть в нескольких последовательных сетках узлов.

Методы решения нелинейных задач в строительной механике приобретают специфические формы. Рассмотрим их на примере нелинейных задач, в которых дискретизация выполняется методом Бубнова-Галёркина в форме МКЭ.

Общим для этих задач является зависимость матрицы жёсткости системы от решения K K(Z). Как было показано в лекции ММ-1, эта зависимость проявляется вследствие того, что коэффициенты дифференциального уравнения задачи несут информацию о свойствах рассчитываемого объекта - его геометрии и механических свойствах материала. В частности, геометрия объекта влияет на вид и коэффициенты уравнений равновесия и уравнений связи деформаций с перемещениями и их производными по координатам, а также в некоторых случаях - на выражения для определения внутренних усилий в сечениях по напряжениям в них. При достаточно больших деформациях проявляется нелинейность их связи с напряжениями. В то же время метод Бубнова-Галёркина предусматривает: а) определение таких перемещений, при которых уравнения равновесия оказываются ортогональными к функциям заданной системы координатных функций; б) выражение усилий через деформации и деформаций через перемещения. Таким образом, указанные уравнения и зависимости оказываются нелинейными. В процессе применения метода дискретизации соответствующие разрешающие уравнения также становятся нелинейными; при этом, как правило, зависимости элементов МЖ от решения записать в виде аналитических выражений невозможно. Эти коэффициенты должны вычисляться каждый раз при новых значениях решения (то-есть, перемещений) по достаточно сложному алгоритму. Такие зависимости естественно называть алгоритмическими.

Уравнение (7.17) принимает вид

K(Z)Z P.

(7.20)

Конкретизируем пути применения изученных нами ранее методов линеаризации (лекции ММ-4 и ММ-5).

1. Метод простых итераций. Запишем (7.20) в виде

Z K(Z) 1 P.

(7.21)

Это по прежнему уравнение, так как оно не разрешено относительно Z , однако, оно приведено к виду, ориентированному на применение метода простых итераций (лекция ММ- 4):

1	(7.22)
Z(i) K(Z(i 1) ) P.

90

В механике деформируемого твёрдого тела, в строительной механике такой итерационный процесс называется методом переменных параметров упругости, так как он предусматривает пересчёт коэффициентов уравнения (7.20) на каждой итерации в соответствии с заданными зависимостями напряжений от деформаций.

Естественно, приведение к виду (7.21) нам понадобилось только для пояснения связи между понятиями метод простых итераций и метод переменных параметров упругости.

При выполнении реальных расчётов записывать решение в таком виде нецелесообразно, нужно просто решать систему уравнений, ставших линейными после подстановки предыдущей итерации в выражение для K(Z)

K(Z(i 1) )Z(i) P.

(7.23)

2. Модифицированный метод Ньютона. Использование обычного метода Ньютона при выполнении реальных расчётов практически невозможно из-за отсутствия аналитических выражений для отмеченных выше нелинейных зависимостей и необходимости их определения в большом количестве точек – узлов системы. Его замена аналогом – многомерным вариантом метода секущих в отмеченных случаях также чрезмерно трудоёмка. В то же время модифицированный метод Ньютона в ряде случаев может быть применён (и действительно применяется).

Запишем (7.20) в виде

F(Z) K(Z)Z P И,

(7.24)

где - нулевой вектор.

Найдём матрицу Якоби левой части (7.24)

F(Z)

(K(Z)Z P) (

K(Z))Z K(Z).

(7.25)

Z

Z

Z

Её значение в точке естественного состояния (Z )

F(Z)

(

K(Z)

) K( ) K( ) K

0

(7.26)

Z

Z

Z

Z

Теперь мы можем записать итерационный процесс модифицированного метода Ньютона

1

(i)

(i 1)

F(Z(i 1) )

(i 1)

(i 1)

1

(i 1)

Z

Z

F(Z

) Z

K

F(Z

)

Z

0

Z(i 1) K0 1(K(Z(i 1) )Z(i 1) P).

(7.27)

Введём обозначения

K(Z) K0

K(Z),

K(Z)Z Pfict , где K0 - матрица

жёсткости системы в естественном состоянии (при нулевых напряжениях), то-есть, на упругой стадии, и подставим их в (7.27):