книги / Проблемы нелинейного деформирования. Метод продолжения решения по параметру в нелинейных задачах механики твердого деформируемого тела

Заметим, что рассмотренные здесь алгоритмы дискретною продол женин в качестве основного элемента содержат систему линейных урав­

нений с расширенной матрицей Якоби ^	Это делает их сходными с ал­
горитмами непрерывного продолжения, и	их алгоритмическая близость

становится еще более очевидной, если системы (1.2.16), (1.2.22), (1.2.26) решать методом ортогоналиэации. Рассмотрим подробное уравнение (1.2.16), а дополнительное условие возьмем в форме (1.2.19);

А*™,)		0.2.28)
< « ( * - ,)/<П1,	Д * <'(; ,)>) = 0.	(1.2.29)
Решение уравнения (1.2.18) можно представить в форме
Д * ^ 1/ =		(1.2.30)
Здесь	произвольный коэффициент,	является реше-
нием однородного уравнения
		(1.2.31)

^- частное решение неоднородного уравнения

Вектор ДАГ° ^ +1^ ортогонален подпространству, натянутому на строки

матрицы У ( Х ^ ку), и может быть найден методом ортогоналиэации, как

и вектор 4X 1	при непрерывном продолжении	(см.	(1.1.28)). При
этом удобно в качестве вектора О принять вектор ДЛГ0^			предыдущего
приближения. Таким образом,
ДХ»«*),> = ог1(7(*™ ), Д * » ” ).			(1.2.33)
Как отмечалось в	§ 1.1, процесс ортогоналиэации	матрицы равносилен

ее представлению в виде произведения левой треугольной матрицы на

ортогональную (1.1.19). В нашем случае			матрица У ( Д Г ^ )	в процессе
ортогоналиэации представляется в виде
г /у (0 т = о (0 у (О				(1.2.34)
(к)>	^ (к) У (кУ
Учитьшая это, помножим систему (1.2.32) на матрицу (Д				) '*, обрат-
ную для Д а)		В результате получим систему уравнений
у (О ду*(»+1) = _ г о (*) I-1 р ( Х ^ 1=			И** .
к (*)а л	(к)	Г '-А (*)'	1 (*>

Так как матрица		У ^	ортогональная, то ее произведение на транспониро­
ванную	(Н (< 0)т	является единичной матрицей. Поэтому непосредст­
венной проверкой легко установить, что				решением уравнения	(1.2.28)
является вектор			'(*) -	(Л)Л	(1.2.36)
1(*)
у *0'+1)_/?/(') чт			(|)
Коэффициент		в (1.2.30) определится		из дополнительного условия
(1.2.29)	в виде
				+>Л .	(1.2.37)
				> )

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

в нашем случае матрица П

является треугольной, что существенно

упрощает ее обращение и делает трудоемкость этой операции равной тру­ доемкости решения только одной системы /и-го порядка. Поэтому исполь­

зование матрицы ( П ^ ) -1 в решении (1.2.36) не существенно увеличи­

вает трудоемкость решения системы (1.2.32) и позволяет при этом ис­ пользовать результат решения системы (1.2.31).

Представление в виде (1.2.30) допускают и все остальные из рассмотрен­ ных выше итерационных процессов. Так, для процесса (1.2.22), (1.2.23) получаем из (1.2.23)

Для процесса (1.2.26), (1.2.27) получаем

д ( ,+ 1) —____г*2		_/■►(*)	у 0 )	д у * ( ' + 1 ) VI /( У(О	д у 0 ( | + 1 ) \
а (*)	— V	^ ( к У	* ( к ) >	^ ( к у а л ( к ) > " ^ ( к У	Д Л (Л) >•

(1.2.39)

Обратим внимание на тот факт, что в случае сходимости всех приве­

денных здесь итерационных процессов, т.е. при	-*• Х (к) , матрица
/ (Л Г ^) сходится к матрице / (Лг(к) ) , которая	играет центральную

роль в процессе непрерывного продолжения решения. Поэтому, сравни­

вая уравнения (1.2.32) и (1.1.8), заключаем, что сходится к

ё Х (к^ й \ . Это позволяет легко комбинировать непрерывное и дискрет­ ное продолжения решения, используя последнее по необходимости для периодического уточнения решения.

Как уже отмечалось, успех применения метода Ньютона — Рафсона во многом зависит от начального приближения. В процессе продолжения

решения по параметру качество начального приближения определяется шагом I по параметру X. Представление решения в виде (1.2.30) позво­

ляет контролировать величину шага I с помощью выражений дли 1

(1.2.37), (1.2.38) или (1.2.39). Для сходящегося итерационного процесса знаменатель этого выражения сходится к скалярному произведению

$(*) ~№Х(к -1)Д/Х, </■*■(*)/<( X).	(1.2.40)
Ясно, что чем меньше шаг, тем ближе друг к другу орты	,)Д/Х

и (1Х(к)1<1\ касательных к кривой К при X = Х*_, и X = X* и гем ближе $(*) к единице. Конечно, 5 (*) может быть вычислено только после того, как итерационный процесс сошелся. Но приближенно судить о значении

$ (Л) можно по знаменателям выражений для <*((/*),)	на первом шаге
(при / = 0):
(</Х{к_ п /4\.	(1.2.41)

Поэтому для того, чтобы шаг / был не слишком мелким и в то же время сходимость итерационного процесса была обеспечена, при практической реализации дискретного продолжения решения обычно достаточно следить за выполненном на первом шаге итерационного процесса условия вида

7„ <(«/*(* _ ,)/</Х, Ц » < ‘ > < у и < 1

(1.2.42)

Из опыта можно рекомендовать для констат уи и уи значения '■'0,8 и ~ 0,9 соответственно. При невыполнении условий (1.2.42) необходимо либо уменьшить, либо увеличить шаг /

1.3. Примеры применения различных форм метода продолжения решения

Здесь мы рассмотрим примеры применения разработанных в § 1.1, 1.2 обобщенных форм метода продолжения решения. Наиболее эффективно эти формы работают, когда множество К решений нелинейной задачи является петлеобразной кривой. Как видно из рис. 1.9, при построении кривой К продолжением по параметру Р мы столкнемся с трудностями

при приближении к'предельной точке В.				Р
Чтобы их преодолеть, необходимо сменить
параметр продолжения. Если в качестве
такого параметра выбрать X и перейти от
Р к X на участке		между	точками А и В,
то точку В мы пройдем без осложнений.
Но вновь столкнемся с трудностями при
приближении к точке С,			которая вместе
с точкой А является предельной при выбо­
ре X в качестве параметра продолжения.
Таким образом, для построения мно­				0
жества	решений	в виде	петлеобразной
кривой	путем продолжения по заданно­

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

Мы здесь пока не касаемся трудностей, возникающих в процессе про­ должения около точки ветвления Г>(рис. 1.9).

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

Уравнение лемнискаты в осях Х х, Х 2 имеет рид
Р(Х) = (Х\ + Х \) 2 —2а21‘М - х Ъ - о ,	(1.3.1)

Эта кривая пересекает ось Х г в точках ± а у /Т и 0. Ниже мы рассмотрим случай а = 1.

Считая Х\	и Х7 функциями параметра X, получаем уравнение продол­
жения (1.1.8)
Здесь / -	расширенная матрица Якоби функции Б ( X ) =Б{Х{,Х 2)-
В развернутом виде уравнение		(1.3.2) принимает вид (произведено сокра­
щение на множитель 4)
Здесь матрица / = [ /ц ./х г ]		представляет собой попросту вектор-строку

Ч

- 1,0

в двумерном векторном пространстве К2. Поэтому в данном частном слу чае нетрудно найти единичный вектор 4Я/4Х, удовлетворяющий уравне нию (1.3.3). Это будет один из векторов

	■Ла	/12
лг	117I	ах “ 117I
</Х	/ч	’ й\
_	1171 _	. 117I

Но мы будем искать решение уравнения (1.3.3) дХ!д\ с помощью общей процедуры ортогоналиэации (1.1.23)

</*/Л = о г1 (7 ,(> ).

(1.3.5)

Интегрировать это дифференциальное уравнение будем, начиная с точки А (рис. 1.10),где

Г = [ч /2 ,0 ]т	(1.3.6)
Начальное значение вектора (? примем в виде
<2= [0, 1].	(1.3.7)

В дальнейшем в качестве () будем принимать предыдущее значение век­ тора ах /ах.

На рис. 1.10 представлены результаты интегрирования задачи Коши (1.3.5), (1.3.6) с помощью явных схем различной точности (т.е. методом непрерывного продолжения) с шагом АХ = 2. При этом рассматривались простой метод Эйлера — результаты показаны точками, модифицирован­ ный метод Эйлера - результаты показаны крестиками и метод Рунге - Кутта 4-го порядка - кружки.

Для того чтобы избежать влияния накопления погрешности в точке ветвления О (в начале координат), эта точка проходилась следующим об­ разом: если начало участка перед следующим шагом продолжения реше­ ния попадало в круг радиуса 0,1 с центром в точке О, то начало участка переносилось в точку, симметричную относительно О. Это позволило сохранить накопившуюся-погрешность.

Накопление погрешности при применении метода Рунге - Кутта с ша­ гом АХ = 2 можно охарактеризовать следующим результатомпосле че­ тырехкратного обхода восьмерки леминскаты кривая пришла в точку с координатами Х\ =1,4112, Х2 =-0,0026. Если бы погрешность не накап­ ливалась, то кривая, конечно, пришла бы в точку А с координатами Хх =

= \/Т « 1,4142, Х2 =0.

Применение алгоритмов дискретного продолжения в § 1.2 при про­ должении решения с тем же шагом АХ = 2 дало практически точное зна­ чение координат точек лемнискаты. Причем затраты машинного време­ ни ЭВМ в этом случае были даже несколько меньшими, чем при интегриро­

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

Пологая арка. Одной из немногих нелинейных задач механики твердого деформируемого тела, допускающих точное решение, является задача о пологой арке. Это делает ее удобной моделью для проверки эффектив­ ности методов решения нелинейных задач. Мы рассмотрим эту задачу в той постановке, в какой она была исследована в работе [107]. Будем также пользоваться, в основном, и обозначениями этой работы, где обоз­ начено: и», и — перемещения оси арки в направлениях координатных осей г , х (рис. 1.11); Л —стрела подъема арки; - интенсивность нагрузки

г,го/

вдоль оси г\ N — продольное усилие в арке, Е - модуль упругости мате­ риала; Р , / - площадь и главный центральный момент инерции попереч­ ного сечения арки относительно оси, нормальной к плоскости арки; / =

=/ (х) —уравнение центральной оси арки до деформации. Введем безразмерные переменные

$ = *//,	И' = *	7=//Л . Ч = Яг 1*1Е1И,
1/ = и/Н,	к = - т 2 /Е1, т = 4ЛРН2		К ' '

В этих обозначениях уравнения, описывающие деформацию арки, концы которой закреплены от перемещений вдоль оси х(ц(0) =«(/) = 0), записы­ ваются в форме

И>1У +кМ"	- к / " ,		(1.3.9)
/ ( / ' и ' + 1	(IV	+ 1 кт = 0, ( . . . ) '= <1 ( . . . Ш	(1.3-10)
Общее решение уравнения (1.3.9) имеет вид
И/(|) = С1со5лД‘? + С2з т ^ Д | + Сз^ + С4 + ^1(1)-			(1.3.11)

Здесь Щ (5) —частное решение уравнения (1.3.9). Если из условий зак­ репления арки найти постоянные интегрирования С1, . . . , С4 и подставить решение (1.3.11) в уравнение (1.3.10), то в итоге получим трансцендентное уравнение, связывающее параметр нагрузки ? с параметром распора арки к в форме

Р(Я.к) =0.

(1.3.12)

Так, для арки синусоидальной формы под равномерной нагрузкой и при шарнирном закреплении концов имеем

7 (|) = —5Ш7Г$,

(1.3.13)

После определения постоянных Сг, . . . ,С* из условий (1.3.14) получаем

^	^со$-\Д %+ (в			51П лД	| - 1^-
1	, 1	+ Т	к		(1.3.16)
- - ( « - * 1)			-----“ *птг|.
2	3	к	- я

Подстановка этого выражения в уравнение (1.3.10) приводит к уравнению, связывающему д и к :

П я. * )“ (Т / ' [ (24 + *)С052 Щ - +

)'Г8",С05 — +

у/к	(1.3.17)
+ со$2 — [я2А(2тг2 - к) + тк(п2 - А)2 ] =
Это уравнение отличается от приведенного в работе [107]	некоторыми

несложными преобразованиями. Его любезно сообщил нам Н.Н. Андриа­ нов. Относительно д уравнение (1.3.17) является квадратным, что поз­ воляет без труда установить зависимость д (Л). Для значений параметра пологости т < 1 эта зависимость образует на плоскости д,к петлеобраз­ ную кривую, которая для т = 0,2 представлена на рис. 1.12.

Мы использовали для построения этой кривой обобщенные алгоритмы непрерывного и дискретного продолжения. Уравнения непрерывного

продолжения представлены в виде
_ (IX	_ Г	ЪР	ЪР	1	_	Г Я1	(1.3.18)
3 —— =	I =	—	, —-	,	* =	.
й \	[	Ъд	Ьк	]		|. к 1

Выражения для ЪР/Ъд и ЪР/Ък несложно получить из уравнения (13.17). Они достаточно громоздки, поэтому мы их здесь не приводим.

	А _______ 1__к__1__а_____________
	>	> "		X
				}
. Л				4/
------ 70
50	С	0 О	100	150 д 200

Рис. 1.12

Явная форма уравнения продолжения получалась решением (1.3.18) методом ортогонализации

<**/</Х = о « ( /, 0 ) .

(1.3.19)

В качестве начальной точкй была взята точка кривой к(д) вблизи начала координат. Само начало координат для такой точки не годилось, так как в нем производные ЪР/Ъц и ЬР/дк содержат деление на к. Начальное зна­ чение вектора О. было принято в виде 0 = [1,0]. Такое задание 0 учиты­ вает, что в начальной точке кривая к (<7) меньше наклонена к оси В даль­ нейшем в процессе продолжения в качестве вектора 0 принимался орт

4Х/4Х для предыдущего значения X.

На рис. 1.12 точками представлен результат интегрирования задачи Коши для уравнения (1.3.19) методом Эйлера с шагом ДХ =4. Крестикам соот­ ветствует результат, полученный модифицированным методом Эйлера с тем же шагом. Кружками обозначены точки, полученные методом Рунге - Кутта 4-го порядка с шагом ДХ = 4 и ДХ = 2. Здесь мы уже имеем практи­ ческое совпадение с точным решением. С целью сохранения накопившейся погрешности переход через точку бифуркации Б осуществлялся также по условию симметрии, как и при построении лемнискаты Бернулли.

Методы дискретного продолжения ( § 1.2) при движении с шагом X = 4 дали результаты, которые совпадали с точным решением с заданной точ­ ностью. При этом машинное время, так же каки е случае лемнискаты Бер­ нулли, оказалось несколько меньшим, чем при использовании метода Рунге - Кутта 4-го порядка.

Трехстержневая ферма. Рассмотрим симметричную деформацию трех­ стержневой фермы, показанной на рис. 1.13, под действием силы Р. Вели­ чинам, относящимся к центральному стержню, мы будем присваивать индекс ”1 ”, а к боковым —индекс ”2”.

Р и с . 1.13

Пусть N1 - сжимающие усилия в- стержнях; Е - модуль Гука при рас­ тяжении; Р ,/( —площадь и моменты инерции поперечных сечений стерж­ ней; I — длина стержней. Оговоримся сразу, что наша основная цель при рассмотрении этой задачи - получить систему нелинейных уравнений, позволяющую продемонстрировать эффективность различных форм ме­ тода продолжения. Поэтому мы не будем стремиться к слишком точному описанию поведения фермы. В частности, мы не будем учитывать измене­ ний углов между стержнями при опускании узла фермы. С учетом этого

уравнение равновесия узла и условие совместной деформации стержней принимают вид

^ 1 + ^ 2 = Л	(13 20)
Д Л =2Д /2.	(1.3.21)
Здесь Д/, - укорочения стержней. Стержни	будем считать неидеальными

и при подсчете их укорочений будем учитывать изгиб стержней под дейст­

вием сжимающих усилий /У,. Пусть х - координата вдоль стержня;		и», (х),
и»,о С*) -	полный и начальный прогибы /-го стержня (/ =1,2).	Тогда с
точностью	до квадратичных слагаемых ((</»у,/</л:)2 = (* >' , ) 2 < 1)	укоро­
чение стержня можно представить в виде

д/' =	+ Т { [(и,,:)2 ~ К Л * *	(ЬЗ-22)
Для описания иэгибных деформаций при действии		сжимающих усилий

/V, неидеального стержня с начальным прогибом н», 0 (х) воспользуемся уравнениями продольно-поперечного изгиба (см., например, [327,329])

/ГУ,ю”+/У,и>, =	;'0,	1=1,2.	(1.3.23)
Будем считать, что начальный проги			ш имеет форму полуволны
синусоиды
н'/о = И;|0 ып тт.х/1.			(1.3.24)
Тогда уравнение (1.3.23) ввиду шарнирного			закрепления концов стерж­
ней (и^'(0) =	=0) имеет решение вида
ттх			(1.3.25)
= IV, 51п ——

Подставив выражения (1.3.24), (1.3.25) в уравнение (1.3.23), получаем связь между амплитудными значениями полного №/ и начального прогибов стержней:

щ - Ц о / о - щ / я р ,

(].3.26)

п2 Е1

/У-

Р’

Здесь ТУ,* — критические значения сжимающих усилий. Из выражения (1.3.22) после подстановки в него прогибов (13.24) и (1.3.25) и интегри­ рования получим

(1.3.27)

д/- б г + 4т (и,? - н/?°)'

Соотношения (1.3.20), (1.3.21), (1.3.26) и (1.3.27) составляют систему шести уравнений относительно шести неизвестных /V,, Д/,, IV,, / - 1 ,2 . содержащую параметр Р. Исключим из этих уравнений укорочения Д/,. / = 1,2, и введем безразмерные переменные

/у,°= /у,//у;, Р° = У У Л /Г Х = и / д

= IV,о//,

а=/У2*/ЛГ: =У2//!, $=РР141и

/= 1 ,2 .

В итоге получаем следующую систему четырех нелинейных уравнений для Ь'? иЛ^° с параметром Р° (здесь и ниже знак ” у безразмерных величин Опущен):

N 1 +сМ2 - Р =О,
- 2аЛГ, + /ЩВ' \ -	2 К?) -	(В'Л - 2 Й^0)] = 0.	(1.3.29)
а д - А Г ^ - й / ю - О ,	^ 2(1 - Ы 2 ) - Щ й =0 .
При Р = 0 эта система имеет решение
Р =О, ЛГ, =N2 =0, Щ = Й/10,		И>2 = 1У20.	(1.3.30)

Отметим также, что если стержни идеальные, т.е. В'ю = И^о =0,то система (1.3.29) имеет четыре точных решения:

	2	Р,	N 2	1	Щ = Щ = 0;		(1.3.31)
\)Г*х = -				= — Р,
2)	^ = 1,		Ыг ={Р-\)1а,			Щ = { 2 Р - Ъ Ш ^ 2 =0;	(1.3.32)
3)	^ ! = ^ - » , ЛГ2 = 1,				Й /,=0, V I =(Р-За)/20;		(1.3.33)
4)^	= 1,	N 2 = 1 ,		Р= 1+ а,		И/? —2 Й/| = (2в —1)/р.	(1.3.34)

Для а = 1, /3 = 100 эти решения и соответствующие им формы деформиро­ вания фермы показаны на рис. 1.14 в пространстве Р. Из этого рисунка видно, что множество решений системы (1.3.29) сложно меняет­ ся в пространстве и имеет три точки ветвления В 1,Д 2| 2?э, причем В 2 и В 3 являются точками вторичного ветвления решений. Такое поведение мно­ жества решений делает невозможным его численное построение с исполь­ зованием единого, параметра продолжения и требует для каждой ветви