Численный метод устойчивого решения задачи Сен-Венана о кручении стержня с произвольной односвязной областью сечения Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2009

Математика и механика

№ 4(8)

Светлой памяти незабвенного Учителя -Павла Парфеньевича Куфарева посвящается

УДК 519.632:531.262

В.В. Соболев

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

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

Ключевые слова: задача Сен-Венана, кручение стержня, гармоническая функция, краевая задача, численный метод, процедура регуляризации, компьютерная программа.

Рассматривается классическая задача Сен-Венана о кручении прямого призматического или цилиндрического стержня с поперечным сечением произвольной формы, скручиваемого моментами силы, приложенными к концам стержня [1-5]. Пусть поперечное сечение однородного по всей длине стержня представляет собой односвязную область В. Влиянием собственного веса стержня пренебрегаем. Поперечные размеры стержня считаются малыми по сравнению с его протяжением в осевом направлении. За ось стержня принимается линия, соединяющая центры тяжести всех поперечных сечений.

Решение задачи Сен-Венана сводится к определению функции кручения Ф = ф(х, у) - гармонической в области В функции, принимающей на границе Г области В значения

Здесь — означает производную в направлении внешней нормали к границе об-

dv

ласти сечения, cos (v, х), cos (v, y) - направляющие косинусы нормали.

Жёсткость при кручении стержня выражается произведением С = цD модуля сдвига ц на модуль кручения (геометрическую жёсткость):

Функция ф определяется граничным условием (1) однозначно с точностью до произвольного аддитивного постоянного, от выбора которого, впрочем, не зависит величина интеграла в формуле (2).

— = ycos (v, х)- xcos (v, y).

dv r

(1)

d

Другая, равносильная форма представления модуля кручения через функцию у = у (х, у), гармонически сопряжённую к функции ф, даётся равенством

Из литературы известны точные и приближённые решения задачи о кручении стержня для самых разнообразных сечений в форме эллипса, правильного треугольника, прямоугольника, кругового сектора, тавра, различных многоугольников и других, сводимых к указанным случаям подходящим конформным отображением [1- 5]. Численные решения краевых задач Дирихле или Неймана, в частности задачи о кручении стержня, для произвольных областей, в том числе клиновидных, могут быть получены методом конформного отображения [2, 5], сведением к интегральным уравнениям [4-7], а также с использованием хорошо известных вариационных методов [7] или получивших своё развитие в последние десятилетия новых методов решения дифференциальных уравнений в частных производных: сеточных, бессеточных, МКЭ, МГЭ, МКГЭ и др. [7-10]. Специфика некоторых из этих методов (конформных отображений, сеточных, МКЭ, МКГЭ) требует учёта конкретного вида области сечения стержня и/или немалого объёма предварительной «ручной» работы для подготовки модели к обсчёту. В других методах затруднены оценки погрешности и/или учёт неустойчивости решения краевой задачи в дискретной постановке, когда граничные значения задаются в отдельных точках границы и притом, как это часто бывает в практических случаях, с некоторой погрешностью. Погрешности могут возникать по разным причинам. В одних случаях это инструментальные погрешности измерения, в других -погрешности интерполяции. Для рассматриваемой задачи кручения граничные значения в (1), (4) не содержат инструментальных погрешностей, однако при интерполяции в межузловых точках границы погрешности возникают.

Поясним сказанное. Для численного моделирования границы Г области В на практике приходится ограничиваться выбором конечного, обычно не слишком большого (из-за ресурсных ограничений) числа т точек (узлов сети) <^1,—, Ст, расположенных на Г или вблизи Г, достаточно полно и точно описывающих Г. При этом в точках границы Г, не совпадающих с узловыми, в качестве граничных значений краевой задачи по необходимости берутся так или иначе интерполированные, приближённые значения.

В простейшем случае кусочно-постоянной интерполяции оценить такую погрешность можно следующим образом. Пусть 81 - положительное число, такое, что граница Г покрывается системой замкнутых кругов радиусом 81 с центрами в узловых точках ^1,.-, Ст , расположенных на Г. В случае выпуклой или полиго-

Функция у , связанная с ф известными условиями Коши - Римана

дф _ ду дф _ ду дх ду’ ду дх ’ должна быть решением краевой задачи Дирихле с граничным условием

(4)

нальной области сечения В за 81 можно взять 81 = — тах 1^,■+1 - С1 (считаем

2 1< ] <т'

Ст+1 = С1). Тогда при замене в точке С еГ точного граничного значения у(0 = |С|2/2 на приближённое у*(0 = |С12/2, где ^^ - ближайшая к ^ узловая точка, возникает погрешность, не превосходящая величины е1 = 81 (Ятах +81 /2). Здесь Ятах - наибольшее из расстояний точек 0,---, Ст от начала координат. Действительно, для ^ еГ имеем

|у(0-у(С ; )| = |а-|С М (|С|+|С ; |)< 2 81 ((тах + Ятах +81 ) = 81. (5)

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

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

1. Постановка краевой задачи

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

замыкании В = В и Г области В, гармоническая в В функция, удовлетворяющая граничному условию (1). Известно, что всякую гармоническую в односвязной области В функцию можно рассматривать как действительную или мнимую часть некоторой регулярной в В функции -^2), 2 = х + 1у . Пусть ф = Яе^. Тогда у = 1т^.

Согласно теореме Уолша [11], регулярную в В функцию ¥, непрерывную в В , можно аппроксимировать многочленом сколь угодно точно равномерно в В : для

любого числа е > 0 существует многочлен

Pn(z) = Z (ak- ibk )zk

k=0

с действительными коэффициентами ak,bk , такой, что max|F(z)-Pn(z)| <е.

zeB

Отсюда имеем max_|у(x,y)-ImPn(x + /y)| <е . Согласно принципу максимума

(x,y )eB

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

max |y(x,y) - ImPn (x + iy)| < е.

(x,y )еГ

Таким образом, функции ф( x,y), у( x,y) можно определить приближённо формулами ф^,у) = фп (x,y) = RePn (x + iy), у(x,y) = yn (x,y) = ImPn (x + iy) с некоторым номером n и коэффициентами ak, bk , выбранными так, чтобы при заданном е > 0 выполнялось неравенство

уп (x,y) -—(x2 + у2)

<є, (x,y) є Г.

(6)

Запишем представления многочленов фп, у п и условие (6) в полярных координатах. Пусть ^1 = |х + іу\ = г, 9 = Л^(х + іу). Имеем

фп (x y) = a0 + Е rk (ak cosk0 + bk sink0);

k=1

уn (x,y) = -b0 +Еrk (bk cosk0 + ak sink0);

k=1

и условие (6) принимает вид

n 1 2

b0 +Еrk (bk cosk0-ak sink0) +— r"

k=1

< Є .

(7)

(8)

(9)

Учитывая, что

получаем

д _ д sin 0 д д . _ д cos 0 д

— = cos 0-------------------------, — = sin 0--------------1---------------

дx дr r д0 дУ дr r д0

- дфп = дуп = Е krk-1 (-bk cos(k -1)0 + ak sin(k -1)0);

дУ ^ k=1

дфп дуп

= Еkrk 1 (ak cos(k-1)0+ bk sin(k-1)0).

(10)

(11)

dx ду k=i

Аналогично предыдущему, представляя функцию ф приближённо гармоническим многочленом вида (7) и учитывая, что

д д . д

— = cos y---+ sin Y— ,

dv dx dy

где y - угол, образованный внешней нормалью к границе Г в точке z е Г с по-

ложительным направлением оси абсцисс, получаем

дф n

—vl = £krk-1 (ak cos ((k -1)0 + y) + bk sin((k -1)0 + y)) , (12)

dv k=1

и граничное условие (1) приводит к приближённому равенству

^krk 1 [akcos((k- 1)0 + Y) + bksin(k-1)0 + y)] =ycosY-xsinY . (13)

k=1

2. Определение функции кручения ф

1) Как известно [12], задача разложения функции в ряд Фурье обладает свойством неустойчивости: при малых вариациях коэффициентов ряда его сумма может измениться сколь угодно сильно. Поэтому даже при незначительных погрешностях, допущенных в определении коэффициентов а0, Ь0,..., ап, Ьп, погрешности в определении фп, уп по формулам (7), (8) при больших п (и, особенно, при г>1), могут стать значительными. В связи с этим для повышения устойчивости решения задачи применим хорошо известную процедуру регуляризации по Тихонову [12].

Задавшись номером п, определим эти коэффициенты в соответствие с условием (13) по методу регуляризующего функционала [12, с. 141], т.е. так, чтобы квадратичная невязка граничного условия (1) по контуру Г была согласована с условием (13) и при этом высокие гармоники многочленов (7), (8) были «подавлены». Именно, потребуем, чтобы было

Ых = Ых(а1,Ь1,...,ап,Ьп) = ст + ХО ^ тш,

1 д \2 где СТ = ст[фп ] =СТ(а1,Ь1,...,ап ,Ьп ) = - {уС°^^ - ^

- квадратичная невязка граничного условия (1); в качестве «стабилизатора» взят функционал

П = ^[фп] = ОЦ,Ь1 ,...,ап,Ьп) = —Г{|э^фп|2 ^ .

ЦГ) {

ЦГ) означает длину границы Г, ds - элемент длины дуги; X, X > 0, - параметр регуляризации, который необходимо определить.

Согласно теории некорректных задач [12], задача на минимум Мх на множестве всех функций фп , удовлетворяющих при заданном 80, 0 < 80 < ст[0], условию ст[фп ] < 80, имеет единственное решение, которое достигается в случае равенства ст[фп ] = 80. Это следует из того, что соответствующее множество функций фп выпукло, а функционал ^[фп ] квадратический.

Записывая необходимые условия минимума дМх дМх

----^ = 0,---^ = 0 (, = 1,2,...,п),

дak дЬ,

приходим к следующей системе линейных алгебраических уравнений (СЛАУ):

Е aC (Х) + Е bD (х) = Fk l=1 l=1

(k = і,..., n). (14)

X alDlk (X) + £blD*kl (X) = Fk 1=1 l=1

Здесь обозначено:

Ckl(X) = Jrk+l-2 ((cos((k- 1)0 + y)cos((l- 1)0 + y)+Xcos(l-k)),

Г

Dkl (X) = Jrk+l-2 ((l cos((k - 1)0 + Y)sin ((/ -1)0 +y)+X sin(l - k)0),

Г

D* (X) = Jrk+l-2 (klsin ((k -1)0 + Y)sin ((/ -1)0 + y)+ X cos(l - k)0),

Г

Fk = k J rk (ycosY - xsinY )cos((k -1)0 + Y)ds,

Г

Fk = kJ rk (ycosY - xsinY)sin((k -1)0 + Y)ds .

Г

Симметричная матрица системы уравнений (14) является матрицей Грама, поэтому она положительно определённая и СЛАУ (14) имеет единственное решение.

Параметр регуляризации X > 0 определим по методу невязки [12] из уравнения

°n () = V (15)

Здесь CTn (X) = a( a1 (х) , b1 (x),..., an (x) , bn (XS) , a1 (X), b1 (X) ,..., an (X), bn (X) -

решение СЛАУ (14) при фиксированном X > 0 .

Величину допустимой невязки 50 > 0 зададим с учётом требуемого уровня

точности е 0 аппроксимации функции дф/ dv на Г: 50 =е;^.

Обычным образом [12, с. 72, 73] доказывается, что an (X) - возрастающая

функция и an (X) ^ с(ш) при X , где

ст(ш) = L1) J(ycosy - xsinY)2 ds.

Следовательно, уравнение (15) имеет (и притом единственное) решение X тогда и только тогда, когда выполняются условия

an (0) <е2 <а(ш). (16)

Если условие an (0) < е2 не выполняется, следует увеличивать n и повторять все вычисления до тех пор, пока не станет an (0) < е2 .

Решение уравнения (15) при выполнении условий (16) можно получить, применяя какой-либо численный итеративный метод, например известный метод хорд.

3. Определение функции у

При известных значениях коэффициентов а1; Ь1,...., ап, Ьп для определения функции у достаточно определить величину Ь0 и воспользоваться формулой (8). Найдём Ь0 из условия минимума квадратичной невязки граничного условия (4):

\2

' ^ ШІП .

Г V к=1 ТІ

Получаем

{( b0 +Еrk ( cosk0-ak sink0) + 2r2j ds ■

2 /7( cosk0-ak sink0) + 2r2 V • (і7)

2L(i ) Г Vk=1 2 j

b0 =-

4. Определение геометрической жёсткости при кручении стержня

Примем приближённо за величину Б значение интеграла

Dn = Я(х2 + у 2 + x “фт - у ~Фп jjdxdy •

в

Применяя формулу Остроградского - Грина, получаем

Бп = ф(-х 2 У - хФп Ух + (У2 х - УФп) ЛУ • (18)

г

Отметим, что величина контурного интеграла (18) не зависит от выбора коэффициента а0 в формуле (7).

Для оценки погрешности при вычислении Б воспользуемся следующей формой представления:

Dn=я(х2 + у2 - х - у dyy^fxdy •

B

Отсюда, учитывая граничные условия (4), приходим к равенству

Dn = 2 jjуndxdy - jj (x2 + y2 )xdy. (19)

B B

Пусть абсолютная погрешность аппроксимации функции у на контуре Г не превышает е : max|у(С)-у„ (С)|<е. Тогда max_|у(x,y)-yn (x,y)|<e и с ис-

СеГ (x,y)eB

пользованием (19) получаем оценку

ADn = |D - Dn\ < 2 jj|y - уn | dxdy < 2eS (B), (20)

B

где S(B) - площадь области B.

Дадим оценку величины е . Пусть Z - произвольно фиксированная точка на Г

и Z j - ближайшая к Z узловая точка. Имеем

|у (О -у n (01<|у(0-у(С j )+|у(С j) -у n (()|+|у n (С j) -y n (C)| <ei + e2 + ез,

где ei = Si (Rmax + Si /2) - ранее установленная оценка (5),

Є2 = max

1< j <m

Уn(ZjHjГ/2 , єз = max|zmax5Jyn (Cj)-yn (z)|

- (апостериорные) оценки, которые при известных коэффициентах ак, Ьк могут быть получены непосредственно с использованием формулы (8). Расчёт величины е2 не вызывает затруднений. Дадим оценку величины е3.

Пусть Д§= Яе (С-С}), Дп = 1т (С-С}), / = (С-С}), г} =|С ;|, 0} С;.

Тогда при |С - Сj | -§1 имеем (0 < 9 < 1)

К (, )-^п К)|=|: V п ( + »«/ )+|; ^П ( +88,е" )

= -9^

S krj1 [acos (k- 1)e j + *)+aksin ((-1)0 j + *)]

k=1

Следовательно, ез < KnS , где

<35і j^-1 (ak| + \bk I)

k=1

Kn = гак Skrk-1 (lakl + lbkI). (1)

1< i<mk=1

Итак, e<e*n, где

e*n = S ((max +Kn + S1 / 2)+e2, (22)

и, следовательно, согласно (20), получаем ADn < Дn, где

Д n = 2e*nS (B). (23)

5. Практическая реализация метода

1. Зададим на границе Г области сечения B стержня произвольно систему точек С1, С2,. ., Сm , образующих на Г сеть узлов так, чтобы замкнутый многоугольник Гт с вершинами в этих точках достаточно точно аппроксимировал область B. Нумеровать точки будем в положительном направлении обхода границы Г (при котором точки области B остаются слева). Будем считать, что угловые точки границы Г, если таковые имеются, включены в данную систему точек. Таким образом, вместо области B везде в дальнейшем берётся близкая к ней (или совпадаю-

щая с ней - в случае полигональной области B) область Bm - внутренность многоугольника (полигон) с вершинами С1, С2,. ., Ст . В качестве длины границы ЦГ) приближённо берётся суммарная длина всех звеньев замкнутой ломаной

С1С 2...С тС1 .

Если исходная область B представляет собой полигон с вершинами A1, A2,..., AP , то в качестве исчерпывающей информации о границе Г достаточно

задать только координаты этих вершин. Однако, если число вершин P невелико, то для повышения точности квадратур необходимо пополнить заданное множество граничных точек дополнительными точками так, чтобы полученная система образовывала достаточно плотную сеть узлов на Гт. Такое пополнение выполняется в автоматическом режиме путём добавления точек, лежащих внутри звеньев ломаной A1 A2...ApA1 . Количество дополнительных точек задаётся оператором по запросу программы. При этом распределение дополнительных точек внутри звеньев ломаной может быть выбрано либо равномерное, либо неравномерное -

со сгущением вблизи вершин А1, А2,..., АР . Будем считать, что пополненное множество образовано вершинами 0, С2, —, Ст . Таким же способом будем обозначать исходное множество граничных точек полигона Вт и в том случае, если пополнение не проводилось.

Если это необходимо (по выбору опции меню), начало системы координат приводится к центру тяжести области Вт и координатные оси - к главным осям инерции. Для этого вычисляются (сведением к контурным интегралам с применением формулы Остроградского - Грина): площадь области Бт, статические, осевые инерциальные и центробежный моменты Ц хёхёу, Ц уёхёу , Ц х2 ёхёу ,

Вт Вт Вт

[[ у2 ёхёу, [[ хуёхёу . По указанным величинам известным образом опреде-

Вт Вт

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

Все криволинейные интегралы по длине дуги на Г вычисляются приближённо как соответствующие определённые интегралы по промежутку 0<5<ДГт). При этом используются значения подынтегральных функций в узлах сети 0 = 51 < s2 <... < $т < sm+1 = ЦГт), соответствующих вершинам ломаной

С1С 2...С тС1. Интегрирование выполняется по специальной, высокоточной квадратуре [13, с. 12], аналогичной известной формуле парабол Симпсона, но не требующей, в отличие от последней, использования обязательно равномерной сети узлов.

2. Задаются положительные числа е0,е4 - допустимые величины абсолютных

невязок граничных условий (1) и (4) соответственно. При заданном номере п формируется расширенная матрица СЛАУ (14) при X = 0 и находится решение СЛАУ. Вычисляются величины сп (0), и проверяются условия (15). Если эти условия выполнены, то совершается итерационный процесс приближённого решения уравнения невязки для параметра регуляризации. Способом «пристрелки» определяется интервал (Х0; Х1) изоляции корня X (0 - X0 < Х1). На каждом к-м шаге итерации методом хорд уточняется значение параметра Xk . С новым значением X = Xк решается СЛАУ и т.д. Условием окончания процесса итераций считается приближённое выполнение равенства (15) с заданной допустимой относительной погрешностью Д0 > 0 : |сп (Хко )/е2 -1 < Д0. Полученные коэффициенты а1, Ь1,..., ап, Ьп используются для вычисления Ь0 (согласно (17) и, далее, значений функции уп в граничных точках ^1, С2,—, Ст . Сравнением с известными граничными значениями у в этих же точках определяется максимальная невязка граничных условий (4). Если эта невязка достаточно мала (меньше заданного е4), то решение задачи приближённого определения гармонических функций ф, у считается законченным. В противном случае номер п увеличивается, в случае необходимости величина е0 уменьшается и все вычисления, начиная с формирования расширенной матрицы СЛАУ при X = 0, повторяются.

3. По найденным коэффициентам Ь0, а1, Ь1,..., ап, Ьп вычисляются апостериорная оценка е*п точности аппроксимации функции у на границе Г и предельная абсолютная погрешность Дп геометрической жёсткости Б согласно (21) - (23).

*

При этом для повышения точности расчёт еп выполняется по более плотному множеству узлов, нежели множество 0, С2, .., Ст . Для этого на каждом звене ломаной , С ]+1 ] (У=1,--,т) берётся дополнительно достаточно большое число (в

нашей программе - до 400) равномерно распределённых точек, и величины Ятах, 81, е2, Кп рассчитываются для новой сети узлов. За окончательное решение принимается найденное путём подбора параметров т, п, е0, е4, X значение Бп с минимальной величиной Дп .

Величина Бп вычисляется, согласно (18), по приближённой формуле

т

Бп =Ё ((( -% Ф*п,; )(+1 -% ) + ( ( -П; Ф*п,; )(( -П; )) .

1=1

Здесь ф*п 1 =(фп (С 1 ) + Фп (С 1+1 ))/2 - усреднённое по двум соседним граничным точкам ^ 1, С1+1 значение многочлена фп ; %, п- прямоугольные координаты точки ^] .

4. Укажем способ определения величины угла у , под которым внешняя нормаль к границе Гт полигона Вт в данной точке С еГт наклонена к оси абсцисс. Для внутренней точки ^ отрезка [^к, Ск+1 ] (к=1,---, т) принимается у = ук , где ук - угол, образованный с осью абсцисс нормалью к отрезку [^к, Ск+1 ]. Для угловой точки 0 границы Гт за величину у принимается среднее значение углов ук-1, ук, найденных для двух смежных отрезков [^к-1, Ск ], [Ск, Ск+1 ] (считаем С0 = С т). Такой выбор угла равносилен определённому сглаживанию границы в точке 0.

6. Опробование метода

Разработан комплекс в составе трёх компьютерных программ [14], реализующих описанный метод. Программа ВоыпёЫеыт подготавливает для дальнейшей обработки необходимые данные: пополняет, если это требуется, множество граничных точек заданного полигона, приводит систему координат к главным осям инерции и рассчитывает граничные значения нормальной производной искомой гармонической функции как решения задачи Неймана. Программа RegulNeum предназначена для решения краевой задачи по методу регуляризации. Полученные в результате работы этой программы коэффициенты гармонических многочленов фп, уп передаются для дальнейшей обработки программой RigidShaftNeum, которая позволяет определить величину геометрической жёсткости стержня при кручении. Кроме того, эта программа может быть использована для построения полей функций ф, у и

касательных напряжений тх =ц05Т/ ду, т^ =-ц05Т / дх в плоскости сечения стержня ( 0 - угол поворота сечений на единицу длины).

Выполнено опробование метода расчётами для различных видов и размеров областей сечений: эллиптическое, треугольное, прямоугольное, в виде кругового сектора, круга с радиальным разрезом, уголка, двутавра и др. Полученные результаты сопоставлены с результатами определения величины D по методу ГЭ (прямая формулировка [8]), получившему в последние годы широкое распространение и признание.

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

1. Модель «Эллипс». Область сечения B - внутренность эллипса с полуосями 5 и 2.

2. Модель «Треугольник». B - правильный треугольник с высотой 3.

3. Модель «Квадрат». B - квадрат со стороной 2.

4. Модель «Шестиугольник». B - правильный 6-угольник со стороной 2.

5. Модель «Круговой сектор-1». B - сектор круга радиуса 1 и раствора 90° (рис. 1).

6. Модель «Круговой сектор-2». B - сектор круга радиуса 1 и раствора 45°.

7. Модель «Круг с щелью». B - круг радиуса 2 с радиальной щелью длиной 1 и шириной 0,1 (рис. 2).

8. Модель «Улитка Паскаля». B - внутренность кривой x = 2 cos t + cos 2t, y = 2 sin t + sin 2t («улитка Паскаля») (рис. 3).

9. Модель «Уголок». B - 6-угольный полигон с 5 прямыми углами; длина большей стороны 3, меньшей - 1 (рис. 4).

10. Модель «Двутавр». B - 12-угольный полигон (рис. 5).

Рис. 3

Рис. 4

В табл. 1, 2 использованы следующие обозначения: Б - точное значение; Бп -приближённое решение; Дп - предельная абсолютная погрешность величины Б;

ЗБп = 1001Б -Бп\/Б - фактическая относительная погрешность (%); Бк - приближённое решение по МГЭ с к линейными элементами.

Т аблица 1

Приближённые значения геометрической жёсткости (Модель «Эллипс», Б = 108,33 [1])

р т Є4 X п Оп Ап Аф Оп 8Оп (%)

64 64 - 0 2 108,73 0,96 0,40 0,37

64 64 - 0 3 108,73 0,96 0,40 0,37

64 128 - 0 6 108,04 0,47 0,29 0,27

64 128 0,00400 0,000117 6 108,07 0,58 0,26 0,24

76 152 - 0 6 108,20 0,43 0,13 0,12

76 152 0,00010 0,000014 12 108,49 0,26 0,16 0,16

100 100 - 0 6 108,50 0,50 0,17 0,15

100 100 0,00010 0,000023 6 108,50 0,50 0,17 0,15

180 180 - 0 3 108,37 0,19 0,04 0,04

180 180 - 0 12 108,37 0,19 0,04 0,04

180 180 0,00001 0,000006 12 108,35 0,16 0,02 0,02

Т аблица 2

Результаты вычисления геометрической жёсткости

Номер модели т п Оп Ап 8Оп (%) О Ок

1 2 3 4 5 6 7 8

1 180 24 108,37 0,15 0,04 108,33 [1] 111,54 (к = 48)

90 180 6 3,1363 0,0226 0,60

2 6 24 3,1222 3,1215 0,0075 0,0045 0,14 0,12 3,1177 [1] 3,1840 (к = 48)

100 8 2,263 0,065 0,44

3 180 10 2,253 0,057 0,00 2,253 [3] 2,283 (к = 48)

180 24 2,253 0,003 0,00

60 10 16,647 0,510 - -

4 192 12 16,576 0,389 - - 16,695 (к = 48)

192 24 16,576 0,102 - -

33 10 0,0828 0,0025 0,36

5 63 18 0,0823 0,0017 0,24 0,0825 [1] 0,1067 (к = 48)

94 26 0,0823 0,0008 0,24

6 195 6 0,0181 0,0000 0,00 0,0181 [1] 0,1222 (к = 48)

Продолжение табл. 2

1 2 3 4 5 6 7 8

7 988 ООО 070 113 24,313 22,582 22,654 - -*) 22,137 (к = 34) 21,894 (к = 60)

8 150 24 53,243 0,435 0,31 53,407 [5] 52,863 (к = 48)

180 30 53,269 0,218 0,26 53,631 (к = 60)

9 192 16 1,587 - - 2,305 (к = 32)

192 28 1,589 - - 1,955 (к = 60)

10 192 16 26 119,89 121,63 - - 127,89 (к = 60)

Примечания: * приведём для сравнения точную величину для сплошного круга: Б = = 25,133; **) неинформативная оценка.

Т аблица 2

Значения функции у и касательных напряжений тгх (модель «Треугольник», Р = 3, т = 90, п = 6, ц = 1, 0 = 1)

№ п/п Координаты точки Приближённые решения Точные решения [9] Относительные погрешности (%)

X У V Т гх V тгх Ау АТх

1 0,0 0,0 0,6668 0,0000 0,6667 0,0000 0,02 -

2 0,5 0,0 0,6459 0,0000 0,6458 0,0000 0,01 -

3 -0,5 0,0 0,6459 0,0000 0,6458 0,0000 0,01 -

4 0,0 0,5 0,6668 -0,5000 0,6667 -0,5000 0,02 0,00

5 0,0 1,0 0,6669 -0,9995 0,6667 -1,0000 0,03 0,05

6 -0,5 0,5 0,6251 -0,7504 0,6250 -0,7500 0,01 0,05

7 -1,0 0,0 0,8336 0,0000 0,8333 0,0000 0,03 -

8 0,5 0,5 0,7085 -0,2498 0,7083 -0,2500 0,02 0,10

9 0,5 1,0 0,8961 -0,5000 0,8958 -0,5000 0,03 0,00

10 -1,5 0,0 1,2290 0,0000 1,2292 0,0000 0,01 -

Как видно из приведённых таблиц, с увеличением числа т узлов на границе, а также с увеличением порядка п гармонического многочлена и уменьшением є4 точность решений, как правило, увеличивается. Для полигональных областей достаточно высокая точность достигается даже без регуляризации (т.е. при X = 0 -по существу, методом наименьших квадратов). Для областей с криволинейными границами, когда дополнительно сгенерированные узлы не лежат точно на границе, применение процедуры регуляризации повышает точность.

Для невыпуклых областей с угловыми граничными точками, при которых внутренние углы больше 180° (см. модели «Круг с щелью», «Уголок», «Двутавр»), расчёт предельных погрешностей Дп приводит к величинам, сопоставимым с величиной Бп («неинформативная оценка»), что несколько снижает ценность метода. Однако и в этих случаях значения Бп получаются близкими к величинам, получаемым по МГЭ.

В сопоставимых условиях, при одинаковых порядках СЛАУ, решаемых нашим методом и МГЭ, величины Бп и Б2п получаются близкими, однако наш метод точнее, особенно для выпуклых областей и областей с гладкими границами. Кроме того, он часто позволяет получать довольно точные результаты даже в тех случаях, когда порядок СЛАУ (14) невелик и при этом намного меньше того порядка

соответствующей СЛАУ в МГЭ, который необходим для получения точности, сравнимой с точностью нашего метода. Высокая точность результатов, полученных в моделях с известными решениями, даёт основание считать метод надёжным. Степень «доверия» к полученному приближённо значению D оценивается по расчётной величине Дп: интервал (Dn -Дn; Dn +Дn) гарантированно покрывает истинное значение D. Метод прост, не требует больших затрат времени для подготовки модели к обсчёту, достаточно эффективен и обеспечивает приемлемую точность решения задачи Сен-Венана. При этом для широкого класса областей указывается гарантированная точность. По затратам машинного времени он сопоставим с МГЭ, несколько уступая последнему, в основном за счёт применения процедуры регуляризации; без её применения быстродействие практически одинаковое: проигрывая МГЭ во времени подсчёта каждого элемента матрицы СЛАУ из-за использования в квадратурах большего числа граничных точек, наш метод выигрывает за счёт симметричности матрицы.

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

ЛИТЕРАТУРА

1. Тимошенко СП, Гудьер Дж. Теория упругости. М.: Наука, 1975. 576 с.

2. Куфарев П.П. К вопросу о кручении и изгибе стержней полигонального сечения // ПММ, 1937. Т.1, вып.1. С.43-76.

3. АрутюнянНХ., АбрамянБ.Л. Кручение упругих тел. М.: Физматгиз, 1963. 688 с.

4. Новацкий В. Теория упругости. М.: Наука, 1975. 872 с.

5. Мусхелишвили Н.И. Некоторые основные задачи математической теории упругости. М.: Наука, 1966. 708 с.

6. Пожарский Д.А. Смешанные задачи теории упругости для составного плоского клина // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2008, №5. С. 36-38.

7. Михлин С.Г. Вариационные методы в математической физике. М.: Наука, 1970. 512 с.

8. Бреббия К., ТеллесЖ.,Вроубел Л. Методы граничных элементов. М.: Мир, 1987. 524 с.

9. Крауч С., Старфилд А. Методы граничных элементов в механике твёрдого тела. М.: Мир, 1987. 328 с.

10. Громадка Т., Лей Ч. Комплексный метод граничных элементов в инженерных задачах. М.: Мир, 1990. 303 с.

11. Walsh J.L. Ueber die Entwickelung einer analytischen Function nach Polynomen // Munchen. Math. Ann. 96, 1926/27. P. 430-436.

12. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1979. 288 с.

13. Соболев В.В., Ищенко Н.В. Численное интегрирование. Методические указания к лабораторной работе с использованием ЭВМ. Ростов н/Д, РГАСХМ. 1999. 28 с.

14. Соболев В.В. Программы численного решения задачи Сен-Венана о кручении стержня произвольного сечения (программный комплекс для ЭВМ). Ростов н/Д, РГАСХМ. За-регистрир. ГОФАП (ВНТИЦ), № 50200802492, 2008. 22 с.

СВЕДЕНИЯ ОБ АВТОРЕ:

СОБОЛЕВ Вадим Владимирович - кандидат физико-математических наук, профессор кафедры «Математика и механика» Ростовской-на-Дону государственной академии сельскохозяйственного машиностроения. E-mail: sobolev@aaanet.ru

Статья принята в печать i8.ii.2GG9 г.