CC BY
44
Вычислительные технологии
Том 11, № 3, 2006
ОБ ОДНОМ ГЕНЕТИЧЕСКОМ АЛГОРИТМЕ И ЕГО ПРИМЕНЕНИИ В ОБРАТНЫХ ЗАДАЧАХ ИДЕНТИФИКАЦИИ УПРУГИХ СРЕД*
И.В. Баранов Донской государственный технический университет, Ростов-на-Дону, Россия
А. О. Ватульян
Ростовский государственный университет, Ростов-на-Дону, Россия
e-mail: [email protected]
А. Н. Соловьев Донской государственный технический университет, Ростовский государственный университет, Ростов-на-Дону, Россия
e-mail: [email protected]
A genetic algorithm holding an intermediate position between deterministic and global probabilistic methods for searching the extremes of a function of several variables is proposed. Numerical code implemented for this algorithm combined with the ACELAN finite element package is aimed to solve the inverse problems of identification of the mechanical properties of isotropic and anisotropic materials. Initial information for the inverse problem can be taken from the displacement field on the free boundary. Results of numerical modeling are presented for determination of various properties of anisotropic material, including density, Young modulus and the Poisson coefficient of isotropic component.
Введение
Задачи об отыскании экстремальных значений функционалов различного вида имеют широкие практические приложения, в частности в оптимальном проектировании [1], в идентификации геометрических характеристик дефектов в твердых телах, в определении физических постоянных, а также и других прямых и обратных проблемах математической физики. Сведение обратных задач к проблемам минимизации целевых функционалов, которые представляют собой невязку между измеренными и рассчитанными характеристиками статического или динамического поведения тела (напряженно-деформированным состоянием, граничными волновыми полями смещений, скоростями распространения упругих волн и т. п.), предполагает наличие надежного и эффективного в смысле быстроты
* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (гранты № 03-07-90411, № 05-01-00734, № 05-01-00690).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
решения средства исследования прямых задач. В качестве такого инструмента могут выступать классические и неклассические варианты метода граничных интегральных уравнений [2, 3] и метод конечных элементов, реализованный в программном продукте [4], который допускает его гибкое использование (в частности, возможность "пользовательского" программирования или совместного использования вычислительных комплексов, как, например, FemLab — MatLab). Как показывает опыт применения таких подходов в задачах реконструкции трещин [5], параметризуемых конечным числом параметров, соответствующие функционалы являются неквадратичными, и поверхности уровня, которые они определяют в области изменения искомых параметров, имеют сложную структуру. В том случае, когда число этих параметров велико либо функционал имеет несколько локальных минимумов или является сильно овражистым, обычные градиентные методы не дают надежного результата.
В последнее время появился совершенно новый перспективный класс вычислительных алгоритмов для решения задач оптимизации — так называемые эволюционные или генетические алгоритмы (ГА) [6-8], использующие для своей работы некоторые принципы развития биологических популяций. Их появление связано с пионерной работой Дж. Холланда [9]. Несмотря на недостаточную теоретическую обоснованность, в задачах многомерной оптимизации эти алгоритмы имеют определенное преимущество перед традиционными градиентными схемами. Они не требуют задания начального приближения, не используют априорной информации о целевой функции и позволяют с большой вероятностью находить глобальные экстремумы в случае ее сложного рельефа. Однако при решении обратных задач идентификации применение этих алгоритмов сдерживается их медленной сходимостью, что в первую очередь, как было отмечено выше, связано с большими вычислительными затратами при определении значений функционала.
Важное приложение обратные задачи идентификации (особенно в режиме "on line") геометрических и упругих свойств конструкций находят в решении проблем неразрушающего контроля. Одной из задач контроля качества изделий является обнаружение внутренних неоднородностей (трещин, включений, полостей). Эта задача может быть разделена на два этапа: на первом определяются геометрические характеристики неоднородности, на втором этапе идентифицируются свойства включения. Проблему первого этапа можно решить методами интроскопии [10]. Исходной информацией для решения задачи второго этапа может служить поле смещений, измеренное на свободной границе тела, при гармоническом или нестационарном воздействии на него. В этом случае задача идентификации может быть сведена к минимизации некоторого функционала, зависящего от параметров включения (плотность, упругие постоянные), в некоторой заданной области их изменения.
Идентификации упругих свойств материалов посвящено огромное количество работ с момента возникновения определяющих соотношений и до настоящего времени. В последние годы к решению этих проблем стали применяться ГА и нейронные сети (NN). В работе [11] находятся упругие постоянные моноклинного кристалла путем минимизации функционала невязки между измеренными и рассчитанными скоростями распространения упругих волн с помощью ГА. Технология нейронных сетей для определения параметров состояния используется в работе [12] при изучении воздействия сферического индентора на тонкие металлические пленки на упругом основании ив [13] при определении эффективных констант в композиционном материале со случайным распределением фаз. В работе [14] с помощью комбинации метода конечных элементов и ГА определяются параметры интерфейсного взаимодействия в композиционном материале, в работе [7] ГА применяются для визуализации дефектов. В работе [15] NN используются для определения акустических
характеристик (коэффициента звукового поглощения, поверхностного импеданса) поли-уретановых пен. В [16] проводится идентификация неоднородных свойств двумерных тел на основе минимизации некоторых функционалов, исходной информацией для этого служит измеренное поле смещений исследуемой области.
В настоящей работе описан оригинальный генетический алгоритм, который сочетает идеи эволюционных алгоритмов с достаточно быстрой сходимостью. Осуществлена его программная реализация, ориентированная на определение упругих характеристик анизотропного тела, механических характеристик локализованных внутренних неоднородностей с использованием кластерной версии конечно-элементного комплекса ACELAN [17].
1. Генетический алгоритм нахождения минимума
функции нескольких переменных в заданной области их изменения
Как правило, задача отыскания минимума функционала сводится к отысканию глобального минимума функции нескольких переменных в некотором конечномерном подпространстве. Последовательность оптимизируемых параметров (особь [7]) записывается в виде хромосомы — конечномерного вектора. Для представления хромосом использовано двоичное кодирование (код Грэя [7], который предпочтительнее стандартного двоичного кодирования в силу сохранения непрерывности бинарной операции) с последовательным расположением генов (векторов) фиксированной длины; при этом каждому параметру соответствует один ген и все гены имеют одинаковую длину. Предполагается, что область изменения (гиперпараллелепипед в Rn, где n — число переменных) каждого из разыскиваемых параметров известна. Область линейно отображается в гиперкуб с ребром длины 2s. Таким образом, каждый параметр переводится в целочисленную точку на отрезке [0...2s], содержащем 2s + 1 точек. Затем это целое число переводится в двоичное представление с помощью кода Грэя.
В алгоритме использованы традиционные генетические операторы, введенные в [7].
1. Двухточечный кроссовер (скрещивание): случайно выбираются две позиции — ¿1 и i2 в хромосоме, хромосомы обмениваются своими частями, которые расположены между ¿1 и i2. В частности, если позиция i2 равна длине хромосомы, то фактически выполняется одноточечный кроссовер.
2. k-точечная мутация. Случайно выбирается позиция в хромосоме, инвертируется значение бита, причем процедура повторяется k раз; в серии вычислительных экспериментов принято k = 3.
3. Инверсия. В хромосоме случайно выбирается позиция ¿, а затем меняются местами части хромосомы, расположенные до и после этой позиции.
Разработанный оригинальный генетический алгоритм состоит из нескольких этапов.
1. Случайным образом (здесь и далее использовано равномерное распределение) генерируется начальная популяция, содержащая N хромосом (нулевое или начальное поколение).
2. Хромосомы сортируются в порядке от лучшей к худшей в соответствии со значением целевой функции на хромосоме. Счетчику i присваивается значение 1.
3. Осуществляется случайный выбор двух хромосом.
4. С вероятностью рс они скрещиваются и порождают двух потомков, причем если скрещивания не происходит, то потомками считаются исходные хромосомы.
5. К каждому из потомков с вероятностью рт применяется оператор мутации.
6. К каждому из потомков п. 5 с вероятностью рг применяется оператор инверсии.
7. Лучшая из этих двух хромосом отбирается и помещается в поколение вместо (N/2 + г)-й хромосомы. Счетчик г увеличивается на единицу.
8. Пункты 3-7 повторяются N/2 раз.
9. Пункты 2-8 повторяются до выполнения критерия остановки.
Критерием остановки всего процесса является стягивание ядра популяции (некоторого процентного числа Q лучших элитных особей) в сферу заданного достаточно малого радиуса Я£.
В предложенном алгоритме в отличие от традиционной рулетки с секторами, площадь которых пропорциональна приспособленности особи, применена "суперэлитная" стратегия, которая обеспечивает, с одной стороны, гарантированное сохранение лучшего из найденных решений, а с другой — достаточно быструю сходимость процесса. Вероятность кроссовера рс выбирается ~ 0.9. Чтобы уравновесить сильное давление отбора, вероятности мутации и инверсии, в отличие от традиционных схем [7], нужно назначать достаточно большими: рт ~ 0.7, а рг ~ 0.2. При этом отметим, что мутация обеспечивает локальный поиск - уточнение, а инверсия — большие броски в области поиска. Таким образом, настроечными параметрами алгоритма являются число особей в популяции — N, число мутирующих бит — к, параметры Q и Я£, а также вероятности рс, рт и рг. По своим свойствам данный алгоритм занимает некое промежуточное положение между детерминированными и глобальными случайными методами поиска экстремума.
2. Постановка обратной задачи идентификации свойств включений в упругих телах
Рассмотрим установившиеся колебания с круговой частотой ш анизотропного неоднородного упругого тела, занимающего конечную односвязную область V с границей Б. Внутри области V имеется ряд непересекающихся подобластей Vq, д = 1, 2,...,N. Упругие свойства в этих подобластях, отличающиеся от свойств основного тела и, возможно, разные для каждой подобласти, обозначим через Уо = V^и^^}. Граница области имеет участок Б0, свободный от напряжений, и может быть описана объединением непересекающихся областей Б = БоиБ\иБ2 (на Б\ заданы кинематические, а на Б2 — силовые граничные условия). На границе внутренних подобластей предполагается непрерывность векторов смещений и напряжений.
Амплитудные характеристики этих колебаний удовлетворяют системе дифференциальных уравнений линейной теории упругости [18]
о,
г],]
р^и2пг
о.
г]
№
Сг]к1Пк,1 г
1, 2, 3, X € Vq,
д
0,1,..., N,
(1)
граничным условиям прямой краевой задачи
Пг1в1 = П( ), = Ог] Ш] ^ = рг, ¿г^о = 0,
(2)
где пг, ц — компоненты векторов смещений и напряжений; сг^ — компоненты тензора упругих постоянных; р(^ — плотность.
В обратной задаче упругие постоянные и плотности в подобластях Уд, д =1, 2,..., N, являются неизвестными, информацией для их определения служит поле смещений, измеренное на границе 50:
игк = и(0)- (3)
На практике измерение поля смещений наиболее просто осуществляется в конечном наборе М точек границы и при К — в значениях частоты колебаний, поэтому соотношение (3) имеет вид
)|^о = М(0)(жт,^г); т = 1, 2, ...М, г = 1, 2,..., К. (4)
3. Сведение задачи идентификации к минимизации функционала
В случае, когда геометрические характеристики подобластей Уд, д =1, 2,..., N, известны, задача идентификации упругих постоянных и плотностей р(<?), д = 1, 2,..., N, сводится к нахождению глобального минимума функционала Е, задающего норму разности (в некотором пространстве Ь) между измеренным (3) и найденным из решения краевой задачи (1), (2) полями смещений, при этом область изменения искомых параметров определяется их физическим смыслом. Пространство Ь определяется исходной дополнительной информацией (свойствами функций и(0)) и методом решения прямой краевой задачи (1), (2). Так, например, в случае условий (4) и решения прямой задачи методом конечных элементов функционал Е может быть представлен в виде
^Ы К 2 А 1/2
ЕЕЕ К0) (Хт ,41р(9)) - ,41Р(9) )к] , (5)
Чт=1 г=1 га=1
при этом область Ш изменения переменных , Р(,?) известна. Отметим, что при известной геометрии включений соотношение (5) определяет функцию многих переменных р(?), д = 1, 2, ...,N, минимум которой может быть найден при помощи описанного выше алгоритма.
4. Примеры численных реализаций
В качестве первого примера программной и численной реализации предлагаемых подходов рассмотрена модельная задача об определении параметров одного упругого изотропного включения в форме круга (центр х0 = 0.07, у0 = 0.04, здесь и далее длины приведены в метрах) в квадратной области ([0.0, 0.1] х [0.0, 0.1]) для плоской задачи теории упругости. Левая боковая сторона квадрата жестко закреплена, к верхней стороне приложено равномерно распределенное нормальное давление (рис. 1).
Находились три параметра: р, Е, V — плотность, модуль Юнга, коэффициент Пуассона соответственно. В качестве материала внешней области выбрана медь. Эксперименты проводились для различных включений — железо, воздух, вакуум.
Решение прямой задачи, моделирующее процесс измерения граничных волновых полей в узлах на правой боковой стороне, проводилось в АСЕЬАМ.
На рис. 2 представлены амплитудно-частотные характеристики (АЧХ) горизонтального смещения средней точки на поверхности измерения, причем цифра 1 соответствует однородной области, а цифра 2 — области со стальным включением. В численном эксперименте выбрана одна "рабочая" частота и = 21.0 кГц (в соотношении (5) К =1) из частотного интервала, в котором наблюдается значительная разница в АЧХ. На рис. 3 на деформированном состоянии области представлены распределения смещений (их слева, иу справа), область включения характеризуется относительной равномерностью цветовой гаммы. На рис. 4 представлены аналогичные распределения в случае, когда область неоднородности представляет собой полость.
Рис. 1.
-— А 1 1 '1
Г 1 П)' Г
1.0
1.0
10
о 0.5
1.5
2.5 СО
Рис. 2.
Рис. 4.
Так как априорная информация касается только геометрии включения, область поиска выбрана в виде Ш = {р,Е, V| [0.0,10.0 ■ 103] х [0.0, 5.0 ■ 1011] х [0.0, 0.5]}, которая включает в себя как "жесткие" и "плотные", так податливые и легкие материалы. В проведенных численных экспериментах выбрано 10 равномерно отстоящих друг от друга узлов измерения смещений на боковой стороне области, свободной от напряжений. Очевидно, что увеличение количества точек измерения (параметр М в соотношении (5)), размера области измерения повысит точность и быстроту сходимости алгоритма, кроме того, выбор этих параметров должен быть согласован с изменяемостью вектора смещений по поверхности измерения, которая зависит от "рабочей частоты". В численном эксперименте ис-
Таблица 1. Результаты численных экспериментов
Параметр Пример
1 2 3 4 5 6
Материал внешней Медь Медь Медь Медь Медь
области
Материал внутренней Железо Вакуум Железо Железо Воздух
области (включения)
Погрешность в геомет- 10-12 0 0 0 0
рии включения, %
Погрешность в измеря- 0 0 0 1 0
емых на поверхности
тела перемещениях, %
Относительный размер 0.1 0.1 0.01 0.01 0.01
включения
Число обращений к 1790 1956 2602 2500 1195
функционалу
Число поколений 100 115 150 150 150
Найдено решение
р (-103 кг/м3) 7.6963 0 7.812 7.858 0.0133 7.86
Е (-1011 Н/м2) 1.9891 0.1443 1.988 1.9992 0.000152 2.0
V 0.2863 0.4999 0.279 0.289 0.4199 0.29
Значение Е 2.560110-6 3.8891 1.2610-11 1.1510-13 1.4810-9
пользовалось 20 особей в одном поколении при вероятностях рс = 0, 9, рт = 0, 7, pi = 0, 2. Результаты экспериментов приведены в табл. 1. Время одного обращения к целевой функции (функционалу) составляло 0.53 с. В колонке с номером 6 приведены точные значения искомых параметров для примеров 1, 3, 4.
Анализ проведенных экспериментов показывает, что предложенный алгоритм за обозримое время позволяет уверенно идентифицировать как достаточно большие (поряд-
ка 10% характерного размера изделия), так и малые (порядка 1%) включения, причем погрешности в определении геометрии включения не слишком сильно влияют на качество получаемого решения.
На рис. 5 представлены лучшие значения р, Е, V (лучшей особи в поколении) и целевой функции Е в зависимости от числа поколений ГА в одном из экспериментов. Анализ результатов, отображенных на этих графиках, позволяет провести дальнейшую модернизацию ГА, внося в него некоторую "интеллектуальную" компоненту как альтернативу случайного поиска при большом номере числа поколений. Это усовершенствование касается сужения области поиска как по длине соответствующих интервалов, так и по количеству искомых переменных.
На рис. 6 и 7 представлены аналогичные (см. рис. 5) зависимости, когда включение представляло собой полость (на рис. 6 — вакуум, число обращений к функционалу (5) равно 1956, на рис. 7 — воздух, число обращений к функционалу (5) равно 1195).
В качестве второго примера рассмотрено определение свойств включения (Ее) в виде эллипса {(х1 — ж10)2/а2 + (ж2 — ж20)2/Ь2 < Л2}, находящегося внутри четверти кругового кольца {г2 < X + х2 < Я2, ж1,ж2 > 0}, состоящего из менее жесткого материала (А1), чем материал включения. Рассматривались установившиеся колебания с частотой f = = 60.83 кГц под действием равномерно распределенного нормального давления
на сторону х1 = 0, при этом сторона ж2 = 0 жестко закреплена. Кроме того, в расчетах принято г = 0.05, Я = 0.1, ж10 = 0.02, ж20 = 0.07, а = 0.01, Ь = 0.005 (все линейные величины даны в метрах). Дополнительной информацией для задачи идентификации характеристик включения р1, Е1, v1 (плотности, модуля Юнга и коэффициента Пуассона соответственно) служило поле смещений, измеренное в десяти равноотстоящих узлах на внешней границе {ж1+ж2 = Я2} кольцевого сектора. В результате 200 поколений с 30 особями в каждом найдены приближенные значения искомых параметров — р = 7.866 -103 кг/м3, Е = 20.01 ■ 1010 Н/м2, V = 0.2902, причем погрешность по сравнению с точными значениями составила менее 1 %, что гораздо меньше той погрешности, с которой эти величины вычисляются на практике. На рис. 8 на недеформированном состоянии области изображено распределение смещения и1, м2, включение выделяется относительно равномерной интенсивностью серого цвета.
В третьем примере определялись компоненты тензора упругих постоянных трансвер-сально-изотропного тела (Ож1 ж2 — плоскость изотропии). В этом случае имеется пять
независимых констант — с11, С12, С13, С33, С55 [18], входящих в обобщенный закон Гука:
Т = СЦ Бj, (6)
где Тг = Бг = м^, % = 1, 2, 3 (по повторяющемуся индексу нет суммирования), Т4 = а23, Т5 = СТ13, Тб = СТ12, Б4 = М2,3 + М3,2, Б5 = «1,3 + «3,1, Бг = «1,2 + М2,1.
В численном эксперименте рассматривались установившиеся планарные колебания закрепленного по левой стороне прямоугольника в рамках плоской деформации под действием равномерно распределенного нормального напряжения на части верхней стороны, узлы измерения амплитуды смещений располагались на верхней и правой сторонах (рис. 9, слева). Частота колебаний f = 5.68 кГц является нерезонансной и расположена между второй и третьей собственными частотами изгибных мод. Векторное поле смещений с распределением его модуля на деформированном состоянии области представлено на рис. 9, справа.
Следует отметить, что в рассмотренной постановке характер напряженно-деформированного состояния не зависит от величины С12, что наглядно иллюстрируется ее поведе-
хЮ
2.8 !
2.6
2.4
10
1
.......„■*****!**■
и '
х 10
13
10
1.6
-
10 20 30 40 50 60 70
10 20 30 40 50 60 70
х 10
-55
10
п......Л-
х 10
'33
0.35
10
0.25
и.
10 20 30 40 50 60 70
110
10 20 30 40 50 60 70
х 10 С
и и—^
1 г 1
10 20 30 40 50 60 70
Рис. 9.
нием в лучшей особи ГА (в отличие от остальных неизвестных, значения которых стабилизируются уже после двадцатого поколения, величина С12 значительно изменяется на протяжении всей эволюции).
Численный эксперимент проводился при следующих параметрах ГА: число особей в популяции — 40, вероятности pc = 0.95, pm = 0.7, = 0.2, число одновременно мутирующих бит — 5, область поиска {[0.1; 4 х 1010] х [0.1; 4 х 1010] х [0.1; 4 х 1010] х [0.1; 4 х 1010] х [0.1; 4 х 1010]}. Анализ сходимости алгоритма указывает на возможность построения ГА с пошаговым сужением области поиска. В табл. 2 представлены результаты поиска погрешности и точные значения искомых констант при 70 поколениях (число обращений к функции = 2653).
Таблица 2. Результаты поиска погрешности
Показатель Константа
С11 С13 С33 С55
Точные значения (-1011), Н/м2 2.627 1.45 2.16 1.29
Найденные величины (■ 1011), Н/м2 2.684 1.546 2.300 1.299
Относительная погрешность, % 2 6.6 6.5 0.7
Разница в погрешностях для найденных констант, по-видимому, определяется выбранной формой колебаний, добиться однородности погрешностей можно использованием в
функционале (5) набора частот, отвечающих различным формам колебаний.
Список литературы
[1] Frey J., Hransbüry J. Saving money through the use of optimization analysis // Water Eng. & Management. 1997. Vol. 144(8). P. 30.
[2] Бревия К., Телес Ж., Вроувел Л. Методы граничных элементов. М.: Мир, 1987. 525 с.
[3] Ватульян А.О., Ковалев О.В., Соловьев А.Н. Новый метод ГИУ в краевых задачах для эллиптических операторов и его численная реализация // Вычисл. технологии. 2002. Т. 7, № 1. С. 54-65.
[4] Васильченко К.Е., Наседкин А.В., Соловьев А.Н. К расчету АЧХ задач об установившихся колебаниях на основе кластерных технологий в ACELAN // Вычисл. технологии. 2005. Т. 10, № 1. С. 10-20.
[5] Баранов И.В., Ватульян А.О. Идентификация трещины в ортотропной среде // Межведомственный сб. науч. тр. Донецк: ИПМ НАН Украины, 2003. Т. 33. С. 133-137.
[6] Goldberg D.E. Genetic Algorithms in Search, Optimization, and Machine Learning. Reading, MA: Addison-Wesley, 1989.
[7] Chiroiü V., Moldoveanu F., Chiroiu C. et al. Application of genetic algorithm in defects visualization // Rev. Roum. Sci. Tech. 1999. Vol. 44, N 2.
[8] Вороновский Г.К. и др. Генетические алгоритмы, искусственные нейронные сети и проблемы виртуальной реальности. Харьков: Основа, 1997. 112 с.
[9] Holland J.H. Adaptation in natural and artificial systems. An introductory analysis with application to biology, control, and artificial intelligence. L.: Bradford Book Edition, 1994. 211 p.
[10] Ковалев А.В. Поисковые технические средства на основе методов интроскопии // Акустические поисковые системы / Специальная техника. 2000. № 2.
http://st.ess.ru/publications/articles/kovalev4/kovalev.htm.
[11] Chiroiu C., Munteanu L., Chiroiu V. et al. A genetic algorithm for determination of the elastic constants of a monoclinic crystal // Inverse Problems. 2000. Vol. 16. P. 121-132.
[12] Hubera N., Tsagrakisb I., Tsakmakis Ch. Determination of constitutive properties of thin metallic films on substrates by spherical indentation using neural networks // Intern. J. of Solids and Structures. 2000. Vol. 37. P. 6499-6516.
[13] Gotlib V.A., Sato T., Beltzer A.I. Neural computations of effective response of random composites // Intern. J. of Solids and Structures. 2000. Vol. 37. P. 4527-4538.
[14] Kang Y.L., Lin X.H., Qin Q.H. Inverse/genetic method and its application in identification of mechanical parameters of interface in composite // Composite Structures. 2004. Vol. 66. P. 449-458.
[15] Gardner G.C., O'Leary M.E., Hansen S., Sun J.Q. Neural networks for prediction of acoustical properties of polyurethane foams // Appl. Acoustics. 2003. Vol. 64. P. 229-242.
[16] Geymonat G., Pagano S. Identification of mechanical properties by displacement field measurement: a variational approach // Meccanica. 2003. Vol. 38. P. 535-545.
[17] БЕлоконь А.В., Бычков А.А., ВлсильчЕнко К.Е. и др. Разработка кластерной версии КЭ комплекса ACELAN // Параллельные вычисления в задачах математической физики: Сб. тр. всерос. научно-техн. конф., 21-25 июня 2004, г. Ростов-на-Дону / Р/на Д.: Изд-во РГУ, 2004. С. 24-35.
[18] Новлцкий В. Теория упругости. М.: Мир, 1975.
Поступила в редакцию 27 января 2006 г., в переработанном виде — 9 марта 2006 г.
