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

62 Вестник СамГУ — Естественнонаучная серия. 2006. №2(42).

МЕХАНИКА

УДК 534.121

ЗАМКНУТОЕ РЕШЕНИЕ ЗАДАЧ ОБ ИЗГИБЕ ЖЕСТКО ЗАКРЕПЛЕННОЙ ПРЯМОУГОЛЬНОЙ ПЛАСТИНЫ1

© 2006 С.А. Лычев^ С.В.Салеев3

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

1. Задача о поперечном изгибе прямоугольной пластины имеет давнюю историю. Первые попытки теоретического исследования изгиба упругой мембраны были предприняты Л. Эйлером, который рассматривал ее как совокупность двух систем струн, натянутых во взаимно перпендикулярных направлениях (1767 г.) Яков Бернулли-младший, используя аналогичные представления, рассматривал пластину как двойной слой брусьев, пересекающихся под прямым углом. Бернулли получил уравнение колебаний, но оно оказалось неверным (1789 г.) [1]. Первые экспериментальные исследования изгиба пластин были осуществлены Хладни (1802 г.) В своих экспериментах Хладни насыпал мелкий песок на стеклянную пластинку и проводил смычком по ее ребру, вызывая колебания. Песок отскакивал от

1 Представлена доктором физико-математических наук профессором Ю.Н. Радаевым.

2Лычев Сергей Александрович (lychev@ssu.samara.ru), кафедра механики сплошных сред Самарского государственного университета, 443011, Россия, г. Самара, ул. Акад. Павлова, 1.

3Салеев Сергей Владимирович, лаборатория высокопроизводительных вычислений Самарского государственного университета, 443011, Россия, г. Самара, ул. Акад. Павлова, 1.

вибрирующих областей и собирался в точках, остававшихся неподвижными, пластинка покрывалась рядом песчаных кривых. Общий рисунок зависел от формы пластины, положения опор и частоты вибрации [2]. Во время своего визита в Париж в 1808 году Хладни продемонстрировал свои опыты перед аудиторией Французской академии наук. Опыты Хладни привели ученых в такое изумление, что они попросили его повторить свои опыты перед Наполеоном. Увиденное произвело на императора впечатление, и он распорядился учредить специальную медаль весом в один килограмм золота и присудить ее тому, кто сумеет дать теоретическое объяснение опытов Хладни [3]. В числе прочих на соискание премии было представлено исследование Софи Жермен. Работа содержала ошибочные промежуточные результаты, но окончательный результат — уравнение изгибных колебаний пластины — оказался верным. Ошибки были исправлены членом экспертной комиссии Ж. Лагранжем, а С. Жермен премирована в 1815 г.

Первое решение задачи об изгибе шарнирно опертой прямоугольной пластины принадлежит Навье, который представил доклад на эту тему во Французскую академию наук в 1820 г. В течение последующих двух десятилетий Коши и Пуассон занимались теорией изгиба пластин исходя из трехмерной теории упругости и полагали, что все определяемые величины могут быть разложены в ряды по степеням малого параметра — толщины пластины. Такие построения приводили к краевым задачам с избыточным числом краевых условий. Г. Кирхгоф в 1850 г. показал, что появление ’’лишних” краевых условий связано с противоречиями в принимаемых гипотезах. Используя вариационный принцип и кинематические гипотезы, теперь носящие его имя, он сформулировал краевые задачи для пластин при различных способах закрепления. Дальнейший прогресс сдерживался отсутствием, за исключением решения Навье, точных решений. В 1899 г. М. Леви предложил разыскивать решение в виде однократного ряда специального вида. Им были получены точные решения для пластин, шарнирно закрепленных по двум противоположным сторонам и произвольно по двум другим. Чуть позже, в 1900 г., Эстанаве показал, что, выполняя однократное суммирование в аналитической форме, возможно преобразовать решение Навье к частному случаю решения Леви [4].

Решение задачи о прогибе жестко защемленной пластины впервые было построено в 1902 г. Б.М. Кояловичем в его докторской диссертации. Краевая задача была сведена к бесконечномерной системе линейных алгебраических уравнений, для которой Б.М. Коялович сформулировал алгоритм приближенного решения и математически обосновал его сходимость [4]. Приближенные методы решения соответствующей прямой вариационной задачи, приводящей ее к бесконечномерной алгебраической системе, были предложены В.Ритцем в 1904 г., а впоследствии развиты Г. Генки (1913 г.) и И.Г. Бубновым (1914 г.) Бубновым были составлены обширные таблицы, сыгравшие значительную роль в инженерном деле. В наиболее общем виде приближенный метод решения задачи об изгибе жестко защемленной

пластины был предложен С.П. Тимошенко на V Международном конгрессе по прикладной механике в Кембридже (1938 г.) Этот метод применим при любом типе нагружения, включая случай сосредоточенной силы [4].

Отметим, что перечисленные методы позволяют получить представление решения в форме разложений по некоторой системе базисных функций, однако коэффициенты этих разложений не могут быть найдены в явном виде и определяются приближенно из решения конечномерного усечения бесконечномерной системы линейных алгебраических уравнений. При этом переход к усечению большей размерности приводит к пересчету всех, в том числе уже найденных коэффициентов [5]. В развитие теории биортогональ-ных систем (В-систем) элементов гильбертова пространства (М.Г. Крейн) Г.Я. Поповым было показано, что с их помощью удается получить явное решение некоторых задач механики (1979 г.) [6]. Решения по-прежнему представляются в рядах, но коэффициенты разложения определяются точно, т.е. в форме конечных выражений. Предлагаемый в настоящей статье метод наиболее близок этим представлениям [6].

2. Упругое тело, толщина к которого мала по сравнению с другими размерами, называют тонкой упругой пластной [4]. В декартовой системе координат {х, у, z} будем рассматривать прямоугольную пластину V со сторонами а, Ь:

Пусть орты I,] лежат в срединной плоскости пластины. Деформирование пластины рассматривается в рамках теории Кирхгофа—Лява [4]: на перемещения и = ш + + кж точек тела наложены связи (кинематические

гипотезы):

где X, ^ — коэффициенты Ламе, после подстановки в него (1) преобразуется следующим образом:

Здесь Vх = V V = д2х + д'Х —оператор Лапласа.

Соотношения (2), (3) определяют специальный вид функционала потен-

V = [0, а] X [0, Ь] X [-к/2, к/Х], к ^ а, к ^ Ь.

и = к - zV ж, ж = ж(х, у), V = I дх + ] ду. (1

С учетом (1) тензор £ малых деформаций принимает вид

а выражение для тензора напряжений Т

Т = Х^£ + XI (I: £),

циальной энергии деформации Х и = ^ Т : £ dхdydz =

V

= ^ [(X + Хц)^2ж)2 - 4ц ((д^ХдуЖ) - (д2^)2) ]zх dхdydz

к Г

2 С7 = — ^ [(X + 2ц)(У2и02 - 4ц ((<92и0(д2и0 - (д2^)2) ] (1х(1у. (4

к/2

= ^ [(Х + Хц)^^)2 - 4ц ((д^Хд^) - (д2^)2)] dхdy J z2dz,

П -к/Х

который после интегрирования по координате z приводится к виду:

/г3

12 I ' 'У/ ' ^ хг* '^у ху'

П

Дальнейшее преобразование функционала потенциальной энергии деформации может быть произведено путем интегрирования по частям правой части (4):

^(д2хум)2 dхdy = £(д2уж)(дхм) dх - £(дх^Хд^) dy + ^(д^Хд^) dхdy.

П дП дП П

На прямоугольной границе дП при условии ж|дП = 0 в зависимости от стороны либо дхж = 0, либо дуж = 0, следовательно,

2 2]

^ [(дхw)(дхw) - (дХуЖ)2] dхdy = 0.

П

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

2 1] = О ^ (y2w)2 dxdy, О = ^ +^2ц)^ ^ ^

П

Располагая (5) и выражением Ж для работы объемных сил К = кк

к/Х

Ж = ^ К  к ж dхdydz = ^ цж dхdy, ц = ^ к dz,

V П -к/Х

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

Ь = и - Ж ЬЬ = 0

П

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

DV2V2w = 0

3. Итак, в прямоугольной области П= [0, а] X [0, Ь] рассматривается краевая задача

V2 Vхw = /, ж|дП =

Здесь / = q/D — заданная, ж — искомая функции, дП — граница области, п — вектор единичной внешней нормали к границе области П.

Полагаем, что искомая функция ж принадлежит пространству £х(П) функций, интегрируемых с квадратом в области П. Будем разыскивать решение задачи (8) в форме ряда, сходящегося по норме £х(П), частичные суммы которого могут быть выписаны в замкнутой форме4. Такое решение подобно классическим решениям Навье и Леви (для шарнирно закрепленных пластин) и по этой причине, в отличие от решения Ритца и Тимошенко, может быть названы замкнутыми.

Краевой задаче (8) соответствует вариационная задача

Ци-’] = ^^У2и^ - fw^

П

где Э обозначает множество допустимых функций, определяемое краевыми условиями (8):

Э с £х

Будем одновременно рассматривать краевую задачу (8) и вариационную задачу (9), учитывая, что функция ж е Э, минимизирующая функционал энергии (9), является решением задачи (8).

Пусть система функций {ж^(х, у)}^ образует счетный базис пространства £х(П). Тогда разыскиваемая функция ж(х, у) может быть представлена в форме разложения

ТО

ж

к=1

Подстановка (11) в (9) приводит к системе уравнений относительно коэффициентов разложения ак:

1=1

П

с/^и,-//щЫ, = 0.

П

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

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

ли базисные функции Wk(х, у) удовлетворяют условию обобщенной ортонормальности

72,,, /„' ,/\тт2.

Wp> = ^ V2Wk(X, у')У2Wp(X, у') йХйу' = Ъkp,

п

то система (12) распадается на двучленные уравнения, а коэффициенты искомого разложения могут быть вычислены следующим образом:

ak =

/г(/, у'

) Wk (X, у') йХ йу'

Из (11), (14) получаем выражения прогибов пластины

/то

К(х, у; х',у')/(х', У') йх!йу', К(х, у; х', у') = ^ wk(x, у^к(х', у'). а к=1

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

4. Для построения Wk может быть использована специальная модификация метода ортогонализации Грама—Шмидта, предложенная Г.Я. Поповым

[6]. Вводятся вспомогательные функции |^к(х, у)}^ е , образующие базис в £2(П). По ним строится система функций ^^(х, у)следующим образом:

1

Wk =

N =

*11 *21 *k1

*1^1 *2k—1   * ?т* . 1

¥1 ¥2 •• ¥k

*11 *21   **1

*lk-l *2k—1    *kk-l

*и *2k   *kk

*пт = ^п, wm>,

N і = 1.

При этом Wk(х, у) также образует базис в £2(П), причем его элементы удовлетворяют условию ортонормальности (13) [6], а выражение для прогибов пластины w(x, у) может быть записано в виде ряда:

w(x, у) = ^

k=l

* *1

\|/1 (х,у) ... Щ(х,у) Г *11

УЛУ\^-1 п ¥1

/(X, у') йХ йу'.

Для реализации вычислительного алгоритма важно, что явные формулы (16) могут быть преобразованы в рекуррентные:

¥1

W1 = , Wk =

У<¥ь Vl)

k-l

4-1/2

<Vk, Щ) - J]<¥k, Wj) Vk - J]<¥k>Wj)w

j=l

k-l

j=l

Таким образом, элементы ортогонализованного базиса записываются в виде

n

wn = ^ ] ftnktyk, k=l

где коэффициенты enk определяются из рекуррентных соотношений

о 1 о bn(h Й21\1/2 R (h 1/2

Pl,= VST Р21 = "|5[Г"?г) • te = p-p7

б. В качестве вспомогательного базиса предлагается использовать систему функций

Ix\ Iy\ / (m + n - l)2 + m - n - l\

где Фп(?) — счетное множество решений двухточечной краевой задачи d4

—ф„ + Л4фй = 0, ф

Собственные значения \n вычисляются как корни трансцендентного уравнения

cosh

Приближенные решения уравнения (20) могут быть найдены численно. Следует отметить, что для значений \n » 1 уравнение (20) асимптотически вырождается в уравнение cos\n = 0, и для значений \n могут быть получены асимптотические выражения:

п

кп = — + ли, и є Z.

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

Функции фп(г), являясь решениями краевой задачи (19), имеют вид:

Фп

+ 8тК(Лп(1 - г)) + со8(Лп) втК^г) - со8К(Лпг) в1п(Хп)],

а коэффициенты bnk могут быть вычислены по явным формулам:

2 с ( b ,4 a 4

Ькіп — J1 (.Рк> Pm)Jі(^7кі Ят) ^кін^2(Рк)^2(Цт) I J рк ^3 qm

Таблица 1

Корни трансцендентного уравнения и их асимптотические

приближения

Номер Корень Асимптотика Абсолютная погрешность

1 4.73004074 4.71238898 -1.76517644Е-2

2 7.85320462 7.85398163 7.77009878Е-4

3 10.99560783 10.99557428 -3.35504373Е-5

4 14.13716549 14.13716694 1.44989660Е-6

5 17.27875965 17.27875959 -6.26556186Е-8

где

Jз(k) при к = т,

1\(к, т) = < 0 при к + т = 2п + 1, п є Z,

J4(k, т) при к Ф т, к + т = 2п, п є Z;

Ык) = е~г}'к ^ш(Хк) - яЩХк))2;

/3(*) = ^е~2Хк1к sec2 16Хк - 64cos(Хk) 8Ї;Ші(Лк) яіп2 Хк+

+ 4Хк ^cos(4Хk) - 2 sec2 Хк - 4«^(2Хк)) + 10 sin(2Хk)-

- 7 віп(4Хк) + cosh(2Хk) (4Хк + 2 вт(2Хк) + sin(4Хk)) |;

Х2Х2 е-Хк-Хт

+ sin Хк !^п Хт (Хк Хк Хт Хт)) ;

Д(к, т) = 8————— ((Хт віп Хш - \к віп \к) втЬ Х„ \к+

Хк4 - Х4т

_ ат _ о.т + ат ^

— — — ш + 2, — ш — — , аш —

17

— + д / — + 2ш

2 V 4

Приведем первые значения коэффициентов /1 (табл. 2) и коэффициентов /2(к) (табл. 3).

Таблица 2

Численные значения коэффициентов /1

к\т 1 2 3 4 5

1 3.18469247 0 -2.47552743 0 -1.93731276

2 11.49464281 0 -4.27889603 0

3 24.72785951 0 -6.08767260

4 42.89628953 0

5 65.99950345

Из последней таблицы видно, что значения /2 для номеров, больших 5, могут быть приняты равными 0.25.

Таблица 3

Численные значения коэффициентов J2

k 1 2 3 4 5

J2(k) 0.25886297 0.24961157 0.25001677 0.24999927 0.25000003

В заключение отметим, что коэффициенты разложения ап (14) для единичной равномерно распределенной нагрузки (f = 1 = const) вычисляются следующим образом:

/ п р

Wn(x', y') dx'dy' = ^ Рпк Yn(x', y') dx'dy' = ab^ |3nkI(pk)I(qk),

к=1 X к=1

ап = J Wn а

1

1

I( Р) = Ф p(t) dt =

а

2e-Xp sin X

X

p

^sin Xp ch Xp - sh Xpj

Численные значения интегралов 1(р), (р = 1,...,9), отличные от нуля, приведены в табл. 4:

Таблица 4

Численные значения коэффициентов I

P 1 3 5 7 9

I(P) 0.42273053 0.18189081 0.11574905 0.084882636 0.067012608

Приведем также формулу для вычисления коэффициентов разложения ап (14) для единичной силы, приложенной в точке X = Хо, y = yo (f = 6(x -

- xo)8(y - yo)):

/П

6(x - xo)6(y - yo)Wn(x, y') dxdy' = ab ^ впк¥k(xo, yo). a k=1

6. В рамках изложенного алгоритма (8), (16), (19) разработана программа в системе символьных вычислений Mathematica 5.2 [9]. Программа предназначена для нахождения аналитических выражений частичных сумм (17), представляющих перемещения пластины. Текст программы в приложении.

Программа определяет новую процедуру DF[p, a, b, Opt], которая возвращает следующие величины, вычисленные для прямоугольной пластины со сторонами a, b (в зависимости от значения Opt):

• матрицу коэффициентов bmn = (ym,yn) (m,n = 1,...,p), если Opt = = Output ^ Determinant;

таблицу коэффициентов ортогонализации |3mn (m = 1,...,p, n =

= 1,..., m), если Opt = Output ^ OrthogonalizationCoefficients;

У

0

У

0

У

0.8

0.6

0.2

Изолинии поля прогибов

1. Равномерно распределенная нагрузка

_§м

Э(75) 1.10

1.05

1.00

0.95

0.90

0.85

Оценка погрешности частичных сумм

1

1. * - —Ху

0 10 20 30 40 50 60 70

Число слагаемых в частичных суммах

Координаты точек: у=0.5,

------х=0.1

х=0.2

------х=0.3

------х=0.4

х=0.5

П

a=b=1

Расчетная схема

0.2 0.4 0.6 0.8 1

2. Сосредоточенная сила, приложенная в центре пластины

Ж

в

1.10

1.05

1.00

0.95

0.90

0.85

Оценка погрешности частичных сумм

і

.-Л _

її >- ~

(

J Г

1

Координаты точек: у=0.5,

------х=0.1

х=0.2

------х=0.3

------х=0.4

-----х=0.5

П

0 10 20 30 40 50 60 70

Число слагаемых в частичных суммах

у////////////////////

f=5

a=b=1

0.2 0.4 0.6 0.8

2. Эксцентрично приложенная сила

1 х Расчетная схема ф=1)

77777777777777777777 _________а=1_______

Ж

вС75)

1.10

1.05

1.00

0.95

0.90

0.85

Оценка погрешности частичных сумм

л V

1 :л /> г-' * и -V. \

1 і

1

Координаты точек: у=0.5,

------х=0.1

х=0.2

------х=0.3

------х=0.4

------х=0.5

П

Число слагаемых в частичных суммах

У////////////////////

0 10 20 30 40 50 60 70

Число слагаемых

f=5

0.2

0.4

0.6

0.8

а=Ь=1

1 Расчетная схема

1

0

1

0

1

0

• последовательность функций wn (п = 1,...,p), удовлетворяющих условию ортогональности (13), если Opt = Output ^ DecomposeFunctions;

• последовательность коэффициентов разложения ап (п = 1,...,p), если Opt = Output ^ DecomposeCoef ficients;

• частичные сумму w порядка п, если Opt = Output ^ Decompose.

Коэффициенты разложения и частичные суммы вычисляются для

• равномерно распределенной нагрузки единичной интенсивности, если Opt = Output ^ Decompose, Loading ^ Pressure;

• сосредоточенной силы, приложенной в точке xo, yo, если Opt = = {Output ^ Decompose, Loading ^ PointForce, PointForceCoord ^

{xo, yo}};

Исследовался изгиб квадратной пластины с единичной стороной и единичной изгибной жесткостью (a = b = D = 1). Рассматривались следующие виды нагружения:

• единичная равномерно распределенная нагрузка (f = 1 = const);

• единичная сосредоточенная сила, приложенная в центре пластины

(f = 8(x - o.5)6(y - o.5));

• единичная сосредоточенная сила, приложенная в точке с координатами x = o.3, y = o.2 (f = 6(x - o.3)6(y - o.2)).

Результаты вычислений приведены на рисунке. Там же приведены численные оценки сходимости частичных сумм для различных точек пластины. Было произведено сопоставление результатов вычисления по предлагаемой методике с расчетами в конечно-элементной среде ANSYS с количеством элементов разбиения Ю4. Наибольшее несовпадение результатов —

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

Литература

[1] Тимошенко, С.П. История науки о сопротивлении материалов / С.П. Тимошенко. М.: ГИТТЛ, 1957. 536 с.

[2] Белл, Дж.Ф. Экспериментальные основы механики деформируемых твердых тел. Ч. 1. Малые деформации / Дж.Ф. Белл. М.: Наука, 1984. 597 с.

[3] Дайан-Дальмедико, Э. Софи Жермен / Э. Дайан-Дальмедико // Scientific American. Издание на русском языке. № 2. 1992 г. С. 60-66.

[4] Тимошенко, С.П. Пластины и оболочки / С.П. Тимошенко, С. Войновс-кий-Кригер. М.:Наука, 1966. 635 с.

[5] Васидзу, К. Вариационные методы в теории упругости и пластичности / К. Васидзу. М.: Мир, 1987. 542 с.

[6] Попов, Г.Я. Биортогональные разложения в задачах механики / Г.Я. Попов // Прикладная математика и механика. Т. 43. Вып. 4. С. 698-708.

[7] Сеницкий, Ю.Э. Изгиб тонкой прямоугольной пластины при различных условиях закрепления на контуре / Ю.Э. Сеницкий // Известия вузов. Строительство. 1998. №6.

[8] Сеницкий, Ю.Э. Об улучшении сходимости спектральных разложений, представляющих решения начально-краевых задач динамики оболочек / Ю.Э. Сеницкий, С.А. Лычев // Численные и аналитические методы расчета конструкций: труды международной конференции. Самара, 1998. C. 189-194.

[9] S. Wolfram. The mathematica book. 5-th ed/S. Wolfram. Wolfram media 2003. 1301 pp.

Поступила в редакцию 17/I/2006; в окончательном варианте — 17/I/2006

CLOSED SOLUTION OF BOUNDARY VALUE PROBLEM FOR RECTANGULAR PLATE FIXED ON BOUNDARY5

© 2006 S.A.Lychev6 S.V. Saleev7

The closed solution of boundary value problem for rectangular plate fixed on boundary is obtained. Typically the solutions of the problem can be obtained by the Rayleigh-Ritz energy method. This method in general gives infinite system of the linear algebraic equations, which is to be solved only numerically. The obtained solution is of the form of expansion based on complete sets of functions satisfying prescribed boundary conditions and special to relation of orthogonality, generated by functional of potential energy of system. The proposed method leads to diagonal infinite system of the linear algebraic equations while the expansion coefficients can be found independently.

Paper received 17/I/2006.

Paper accepted 17/I/2006.

5Communicated by Dr. Sci. (Phys. & Math.) Prof. Y.N. Radayev.

6Lychev Sergey Alexandrovitch (lychev@ssu.samara.ru), Dept. of Continuum mechanics, Samara State University, Samara, 443011, Russia.

7 Saleev Sergey Vladimirovich, Laboratory of High Performance Computing, Samara State University, Samara, 443011, Russia.

ПРИЛОЖЕНИЕ

A [n_] : = A /. FindRoot [Cosh [A] Cos [A] == 1, {A, — + n nj ,

WorkingPrecision^ 100, PrecisionGoal^ 95, MaxIterations ^ 200J

F [x_, A_] : = -A (Sin [A (1 - x)] + Cosh [A] Sin[A x] -

Cos [A x] Sinh [A] + Sinh [A

DF [n_, a_: 1, b_: 1,

Opt_: {Output ^ Decompose, Loading ^ Pressure, PointForceCoord ^ {1 / 2, 1 / 2}}] : = Module [{L, P, Q, s, k, m, J, J1, J2, J3, J4, D1, D2, T1, T2},

L = Table [A [k], {k, 1, n}] ;

{P, Q} = Transpose[

r r 1 7 _ f s + 3s s2 + s.. ..

Table [s = Floor [----+ -A---+ 2 k J ; {------ - k + 2, k - ----j, {k, 1, n}JJ;

1

J = Table [if [k == m, -- -2 A A Cosh [A] 2 Sin [A] (4 Cosh [2 A] Cos [A] 3 +

L L 16

3 Cos [A] + 16 A Sin[A]3 - 7 Cos [3 A] - 16 Sin [2 A] Sinh [A]) / . A ^ L[[k]],

r A2 m2 -A-M

If[oddQ[k + m] , 0, 8 - ((m Sin[M] -A Sin[A]) Sinh[M] Sin[A] Sinh[A] +

A4 - m4

Sin[A] Sin[M] (A Sin[A] Sinh[A] -m Sin[M] Sinh[M])) /.

{A ^ L[[k]] , m ^ L[[m]]}JJ , {k, 1, n}, {m, 1, k}J ;

J1 = Table [If[m > k, J[[m, k]] , J[[k, m]]] , {k, 1, n}, {m, 1, n}] ;

J2 = Table [-2 A (Sin [A] - Sinh [A]) 2 / . A ^ L[[k]] , {k, 1, n}];

J3 = Table[

2

---- J1[[P[[k]], P [[m]]]] J1[[Q[[k]], Q[[m]]]] + KroneckerDelta[k, m] J2[[P[[k]]]]

a b

J2[[Q[[m]]]] (-b- L[[P[[k]]]]4 + -.a.- L[[Q[[m]]]]4|, {k, 1, n}, {m, 1, n}J;

If [(Output /. Opt) = = = Determinant, Return [J3]];

D1 = VTable[Det [Take[J3, {1, , {1, k} ] ] , {k, 1, in}] ;

D2 = Prepend [Table [1 /

Table [( -1) k+m Det [Drop [Take [J3, {1, k - 1}, {1, k}], {}, {m}]] , {m, 1, k}], {k,

2, n}], {----1----jJ;

1 D1[[1]] Ji

If[(Output /. Opt) = = = OrthogonalizationCoefficients, Return[D2]];

T1 = Table [ V" D2 [[m, k]] f[ -- , L[[P[[m]]]]J f[ , L [[Q [[m]]]]J , {m, 1, n}J;

k:! a b

If[(Output /. Opt) = = = DecomposeFunctions, Return[T1]];

If [(Loading /. Opt) = = = Pressure,

_ 2 -A Sin [A] (Sinh [A] - Tan [A]) _

J4 = Table [---------------------------------/. A ^ L[[k]] , {k, 1, n}J ;

A

m

T2 = Table [a b ^ D2 [[m, k]] J4 [[P[[k]]]] J4 [[Q [[k]]]] , {m, 1, n}JJ ;

k=1

If [(Loading / . Opt) = = = PointForce, J4 = PointForceCoord / . Opt;

T2 = Table [T1 [[m]] / . {x ^ J4[[1]], y -> J4[[2]]}, {m, 1, n}]] ;

If[(Output /. Opt) = = = DecomposeCoefficients, Return[T2]];

n

If [(Output / . Opt) = = = Decompose, Return [^^ T1 [[k]] T2 [[k]]JJ