© Ф.Пинежанинов
К списку статей
Зеркало: http://pinega.da.ru/
Для использования метода конечных элементов вообще-то не нужна
теория оболочек и стержней. По существу нужна технология построения матриц
жесткости, масс и демпфирования для конечных элементов как геометрических
образов на многообразиях в трехмерном пространстве и тогда можно будет
рассчитать оболочки.
Под многообразием, как обычно в математике, будем понимать
обобщение понятия линий и поверхностей без особых точек. В математической
энциклопедии отмечается, что это понятие может быть использовано во всех
ситуациях, когда рассматриваемые объекты могут быть представлены как
отображения, зависящие от параметров. Точками многообразия могут быть объекты
любой природы - точки в геометрическом смысле, матрицы, состояния механических
систем и так далее.
Конечный элемент с двумя параметрами всегда может быть
использован для трехмерного описания поверхности как объекта без толщины, а
наличие гладкости обеспечивает однозначное существование нормали и,
следовательно, если каждой нормали сопоставим пару чисел, то можем наделить эту
поверхность толщиной, а значит, это многообразие может порождать объем. Примерно
так, как используется наследование в объектно-ориентированном
программировании.
В объеме может действовать нормальная механика сплошной среды и
к ней применяться трехмерный метод конечных элементов. Построение матриц
жесткости для трехмерных элементов вполне строгий процесс. Другими словами,
поверхностные элементы вполне могут использоваться только как удобный способ
описания, а практические вычисления можно осуществлять на трехмерных задачах,
используя перемещения в качестве неизвестных параметров. Для гладких
поверхностей достаточно часто так и делается. Как замечает Сьярле " в инженерной
практике этот метод в ряде случаев конкурирует с методами, непосредственно
основанными на двумерной постановке".
Но есть и проблемы у этого подхода. При негладких поверхностях
непросто обеспечить стыковку на линиях излома, в слоистых композитных материалах
никаких ресурсов и терпения не хватит, если моделировать каждый слой, непросто
стыковать, например, две слоистые оболочки со слоями разной толщины. Примеры
можно долго продолжать. В стержнях естественно те же проблемы.
С другой стороны в инженерных теориях оболочек, основанных не
только на перемещениях, но и углах поворотов эти проблемы решаются очень изящно.
Один недостаток - обилие теорий и гипотез, сложный математический
инструментарий, основанный на криволинейных координатах, не позволяет быть
уверенным в надежности получаемых результатов. Оболочка уже тонкая или еще
толстая, например. Да и в самих оболочечных конечных элементах предложений
много, а это говорит о том, что еще не найдено приемлемого решения проблемы. В
расчетной практике, даже Джон Свенсон, основатель Ansys, призывает к
осторожности при использовании квадратичных оболочечных элементов.
Множество
подходов создает "конструктивную неопределенность" как сказали бы инженеры, или
"диалектическое противоречие" ответили бы философы.
Критиковать то легко, "каждый мнит себя героем, видя бой со
стороны" писал Шота Руставели, поэтому, вернемся к многообразиям и голой
абстракции, потому, что только они позволят понять существо проблем, через
взгляд издалека и отсечение деталей, которые не позволяют охватить всю
совокупность явления.
"Ибо не добыча, а само понятие является самым
главным...", так вроде говорили герои фильма о национальной охоте.
Как отмечает Работнов - наиболее надежный способ построения
приближенных теорий состоит в том, что за основу берутся вариационные принципы и
делаются некоторые предположения о характере распределения перемещений, потом на
этой основе получаются дифференциальные уравнения и необходимые для их
решения граничные условия.
Принцип естественно не вызывает возражений,
но в методе конечных элементов не нужны уравнения, а следовательно и условия,
ими порождаемые. Нам нужны матрицы и возможность их обращения, а, следовательно,
условий достаточно таких и столько, чтобы матрицы были не вырождены.
Два параметра порождают поверхностный элемент. Доопределив его
парой чисел, породим объем или оболочку, на которой из-за деформации возникнет
поле перемещений, используя перемещения в качестве параметров, вариационные
принципы и соотношения трехмерной теории упругости получим матрицы жесткости,
как линейные операторы, действующие на поле перемещений.
С другой стороны, поверхностный элемент в процессе деформации
породит многообразие перемещений точек поверхности и углов поворота, которые
можно использовать в качестве параметров. В сочетании с инженерными теориями
оболочек, основанными на гипотезах поведения и вариационными принципами,
построим матрицы жесткости поверхностных элементов, как линейные операторы,
действующие на перемещения и углы поворота поверхности.
Применяя обычную конечно-элементную технику, можем решить обе
задачи и в силу линейности задачи теории упругости это должно быть одно решение.
Вспоминая, что по Гильберту существование в математике - это непротиворечивость
и если, выходя из одной точки и идя разными путями, приходим снова в одну точку,
значит, получаем замкнутую область. Ее можно перерезать тропинками, соединяющими
одну линию с другой и, встречая сложности на одной линии, переходить на другую,
в результате все равно придем куда надо, если будем "гнуть свою линию", как
советовали в чудесном фильме о брате.
Примерно так рассуждал автор, гуляя по тропинкам садика и
удивляясь причудливому сочетанию ограды Зимнего дворца, бывшему ресторану, в
котором поп Гапон готовил шествие к царю, памятника Кирову, домику, где было
принято решение об Октябрьской революции и парку в который "ездили гулять
Обломов со своим энергичным другом".
"Идея должна быть в корпусе, а не в принципе" скажут
современные пост модернисты и будут правы. Поэтому перейдем к делу и установим
связь между полем матриц жесткости трехмерных оболочек и стержней и
многообразием матриц жесткости поверхностных и линейных объектов. Заодно и
матриц масс и матриц демпфирования.
Матрицы для трехмерных
объектов
Конечные элементы обычно используются для построения
характеристических матриц в задачах механики деформируемого тела и
математической физики.
Рассмотрим технологию их использования на примере матрицы
жесткости в теории упругости. Эта матрица обычно строится на основе функционала
упругой энергии для конечного элемента.
Для ее практического вычисления надо осуществлять
интегрирование по области конечного элемента. Область обычно имеет произвольный
вид, поэтому используют параметрическое отображение на область порождающего
элемента и интегрирование осуществляют уже на нем. Делается это обычно так:
Поскольку деформации определяются через производные, то
необходимы соотношения переводящие производные в глобальной декартовой системе
координат в локальную систему координат конечного элемента. Они строятся обычным
образом на основе отображения области порождающего конечного элемента на область
реального конечного элемента.
Рассматривая компоненты как сложные функции можно записать
связь между дифференциалами в этих системах координат:
Квадратная матрица,
осуществляющая линейную связь между дифференциалами в различных системах
координат, называется матрицей Якоби отображения или производной Фреше
отображения и обычно обозначается J или, иногда, чтобы подчеркнуть ее структуру
построения - градиентом отображения
 .
Из математического анализа известно, что ее определитель,
называемый якобианом, дает искажение малых объемов при замене переменных под
знаком интеграла, а так же знак якобиана говорит об ориентации бесконечно малого
объема. Как содержательно использовать эту информацию при практическом
интегрировании, подробно описано в статье "Интегрирование конечных
элементов".
Обратим внимание на геометрический смысл матрицы Якоби.
Зафиксируем, например, h,z и будем считать переменной только
x. Тогда объем
выродится в линию, а матрица Якоби выродится в вектор, который будет касательным
к параметрической кривой .
Аналогичная ситуация и в двухпараметрическом случае по
столбцам матрицы Якоби будут векторы, касательные к поверхности описанной двумя
параметрами . Соответственно в трехмерном случае по столбцам стоят касательные
векторы в глобальном декартовом базисе к координатным плоскостям пространства
порождающего конечного элемента.
Деформации определяются через частные производные от вектора
перемещений и если его компоненты рассматривать как сложные функции
то можно записать связь между производными вектора перемещений
в локальной системе координат и в глобальной декартовой системе в виде
Или в символической записи, в которой вектор всегда понимаем
как столбец:
Здесь  -
транспонированный градиент перемещений, или транспонированная производная Фреше
перемещений. Транспонируем это матричное равенство, получим обычное правило
цепочки дифференцирования сложных отображений . Если матрица Якоби не вырождена,
то можно записать и обратное преобразование через матрицу обратную к матрице
Якоби .
Это ключевое для метода конечных элементов соотношение,
позволяющее записать тензор малых деформаций через градиент деформаций в
символическом виде и позволяющее практически вычислять интегралы
через вычисление этих матриц в точках интегрирования численных методов.
Вычисление интегралов для трехмерных элементов практически
стандартная процедура и изложена в любой книге по конечным элементам.
Матрицы для
двухпараметрических многообразий - оболочек
В природе и технике достаточно часто встречаются объекты, в
которых существует одно направление, в котором размер значительно меньше, чем в
ортогональных ему направлениях. Такие объекты, которые имеют две
поверхности, расстояние между которыми или постоянно или слабо изменяется, и
называются оболочками. Эти объекты удобно описывать с помощью некоторой
ориентированной поверхности и толщины, то есть размера в направлении
ортогональном поверхности.
Если поверхность проходит точно по середине между
поверхностями, ограничивающими оболочку, то она обычно называется срединной
поверхностью. Это не обязательно и годится любая близкая поверхность, называемая
опорной или несущей, и пара чисел, которые задают начало и конец оболочки на
нормальном к опорной поверхности векторе, исходящем из опорной поверхности и
направленном в положительном направлении относительно ориентации
поверхности.
Аналогично и в конечных элементах. Двухпараметрические конечные
элементы с помощью преобразования
могут порождать поверхностные элементы в трехмерном
пространстве. В таких элементах всегда неявно подразумевается и третье -
нормальное к поверхности направление, положение объекта вдоль которого можно
определить парой чисел. Это может быть верхняя и нижняя граница поверхностей
оболочки или положение середины и толщина.
Наличие такого особого направления позволяет принимать
интуитивно обоснованные предположения о характере распределения перемещений и
(или) напряжений вдоль нормали к опорной поверхности.
Например, если возьмем пластинку жевательной резинки и на
противоположных сторонах пластинки нанесем продольные царапины одинаковой длины,
а потом согнем пластинку, то увидим, что длина одной царапины увеличится, а
другой уменьшится. Если изгиб осуществляется относительно срединной поверхности
и упруго, то понятно, что перемещения должны быть симметричными и иметь разные
знаки, то есть функция описания перемещений вдоль толщины должна быть
кососимметричной функцией.
С другой стороны если просто будем растягивать пластинку, то
естественно удлинения царапин будут одинаковые, то есть функция, описывающая
перемещения вдоль толщины должна быть симметричной. Вывод таков - функция может
быть произвольной, но, по крайней мере, достаточно плавно меняющейся или точнее
гладкой. Другими словами несколько членов ряда Тейлора должны хорошо описывать
изменение перемещений вдоль толщины, относительно срединной точки.
Пусть оболочка деформировалась, для опорной поверхности это
приведет к тому, что для произвольной точки оболочки изменится положение и
направление нормали.
Разницу между деформированным и начальным состоянием точки
опорной поверхности удобно описать с помощью вектора перемещения u и вектора
углов поворота q.
В теле оболочки на линии совпадающей с направлением нормали
возникнут перемещения, которые можно связать с движением опорной поверхности
следующим образом:
Здесь  - вектор перемещения произвольной точки
оболочки,  - вектор перемещения согласованной с ней точки опорной
поверхности,  - вектор углов поворота опорной поверхности,  -
единичный вектор нормали,  - ориентированное расстояние от точки опорной
поверхности до средней точки оболочки, z - ориентированное расстояние от
точки опорной поверхности до произвольной точки оболочки. Другими словами
движение произвольной точки оболочки можно разложить на три движения: движение
согласованной точки поверхности, движение срединной точки оболочки относительно
точки поверхности и вращение произвольной точки оболочки, относительно срединной
точки. При этом понимается, что все точки лежат на нормали к
поверхности.
Как отмечалось ранее  - должна быть
кососимметричной функцией.
Простейший случай, когда используется линейная функция  ,
приводит к гипотезе плоских сечений, но можно принять из каких-то других
соображений или экспериментов и другие гипотезы о кососимметричной функции
f.
Пусть в некоторых точках оболочки лежащих на прямой линии с
нормалью n определены перемещения, тогда для каждой точки с номером i
можно записать множество соотношений  .
Если этих векторных
соотношений не меньше 2, то можно найти два неизвестных вектора ,
линейно связанных с перемещениями оболочки, при этом если известна точка, то
вектор n тоже определен однозначно из соотношения  , где столбцы
матрицы векторы касательные к опорной поверхности. Поэтому, если обозначить их
за  можно вычислить
 нормальный к поверхности вектор единичной длины . Детали
построения матрицы аналогичны трехмерному случаю. Для удобства манипуляций
введем обозначение
 . Небольшой код дает для одной точки оболочки
три линейных соотношения относительно шести неизвестных {u,v,w, Qx,Qy,Qz}:
 =>
Представим эти соотношения в матричном виде удобном для
составления системы линейных уравнений:
 => .
Таким образом, устанавливается линейная связь между
перемещениями всех точек реальной оболочки лежащих на линии с нормалью n
с перемещением и углами поворота одной точки пересечения этой линии с опорной
поверхностью. То есть для каждой точки опорной поверхности имеем систему
линейных уравнений  , где i номер точки принадлежащей оболочке.
После такой подготовки перейдем к конечным элементам. Пусть
имеем описание геометрии оболочечного элемента в виде ориентированного
нумерацией поверхностного элемента, то есть координат узлов элемента опорной
поверхности и с каждым узлом ассоциированы два числа, определяющие толщину
оболочки.
На основании этой информации построим для каждого узла
поверхностного элемента нормали и затем трехмерный элемент при этом грани
трехмерного элемента должны быть подобны поверхностному элементу, а ребра на
нормалях могут иметь произвольное число узлов, но не меньшее двух. Для этого
трехмерного конечного элемента построим матрицу жесткости, как описано ранее по
стандартной процедуре и этот этап обычно не вызывает трудностей. Как действовать
дальше, рассмотрим на примере матрицы жесткости.
Проблему минимизации квадратичного функционала можно
свести к системе линейных уравнений, которая получается из сложения отдельных
систем для каждого конечного элемента C U = F, далее будем считать, что эта
элементная система построена для трехмерного элемента. Здесь вектор F может быть
основан на массовых силах и (или) начальных деформациях и температурных
воздействиях. Объединение элементов в систему и учет различных условий
подробно описаны в статье "Стационарные точки и условия на переменные".
Кроме того, есть возможность точнее учесть распределенные
моменты, возникающие в оболочке от действия поверхностных сил. Как показано
ранее, справедливы соотношения связывающие перемещения трехмерного элемента и
перемещения и углы поворота для поверхностного элемента U = A Q. Подставим их в
C U = F и получим C A Q = F. После этого умножим слева на
транспонированную матрицу A и получим  .
По существу, произошел переход к новым обобщенным перемещениям
или новым параметрам и обобщенным силам, как это обычно делается в
аналитической механике.
Теперь квадратную и симметричную, если C была симметричной,
матрицу  можем рассматривать как матрицу жесткости оболочечного элемента, а
вектор  как элементные силы, приведенные к поверхностному элементу. Следуя
примеру физиков, соединяющих пространство и время в одних операторах, можно
соединить геометрию и силовые воздействия в рамках одной матрицы, как это
делалось в статье "Стационарные точки и условия на переменные".
Внимательный читатель подумает, и двух точек на нормали
хватило бы для преобразований, "все то зачем" как сказали в популярном фильме о
русской охоте.
"Дело в следующем" если хотим, чтобы была возможность реального
учета нелинейности перемещений по толщине  , и соблюдались
условия кососимметричности, надо не менее четырех точек. С другой стороны мы
можем рассчитывать, например вафельный торт и нам нужна возможность построения
столбика из трехмерных элементов с различными свойствами, возможно в слоях еще
будут орешки, которые учтем, как описано в статье "Осреднение свойств в
конечном элементе" и все это надо смоделировать одним оболочечным элементом с
эквивалентными жесткостными свойствами. Для этого и пригодится возможность
учесть и усреднить в смысле наименьших квадратов свойства приведением к
оболочечному элементу. Понятно, что нет препятствий для применения конденсации
внутренних узлов при использовании техники типа суперэлементной при условии
обеспечения согласованных с оболочечным элементом границ. Такой слоеный пирог из
трехмерных элементов не обязательно должен быть сцеплен по всей поверхности, что
может пригодиться в расчетах при расслаивании слоистых композитов или при
возможности проскальзывания арматуры в наполнителе.
Есть еще одна интересная возможность для соотношения  . При получении в результате решения общей задачи на
оболочечных элементах информации о перемещениях можем снова вернуться к
трехмерному "пирогу" и, учитывая малую толщину и принцип Сен-Венана о
локализации особенностей, можем в середине элемента уточнить гипотезу о
распределении толщины и снова решить задачу. То есть существует итерационный
процесс, позволяющий из первоначальной простейшей гипотезы найти ассимптотику к
истинному распределению перемещений поперек оболочки.
С другой стороны, всегда возможно использовать простейшую
кососимметричную линейную функцию и за счет корректирующего коэффициента,
полученного на основе экспериментов, или известных аналитических решений
произвести подстройку матрицы жесткости оболочки. Несложно этот коэффициент
сделать функцией отношения толщины и характерного размера, или, что практически
то же самое, на основе диагональных членов матрицы Якоби, и таким образом
обеспечить "стирание граней" между тонкими и толстыми оболочками.
То есть, нам представляется возможность смешивать разные
переменные в рамках одного элемента, в части его трехмерные свойства, а переход
его к обобщенным переменным осуществлять только на некоторых группах узлов. Это
позволит легко и плавно соединять трехмерные элементы с поверхностными
элементами, через такие элементы: с одного края оболочки, а с другого
трехмерные.
Матрицы
для однопараметрических многообразий - стержней
Под стержнями обычно понимается трехмерное тело, поперечный
размер которого мал по сравнению с длиной линий образующих боковую
поверхность.
Построение конечных элементов для криволинейных стержней
естественно проводить аналогично оболочечным элементам, ассоциируя со стержнем
опорную линию и поперечное сечение.
При этом стержневой конечный элемент в трехмерном пространстве
можно рассматривать как однопараметрическое отображение порождающего одномерного
конечного элемента:
Как
обычно выражение для полного дифференциала дает нам вектор касательный к опорной
линии, который определяет плоскость ортогональную опорной линии:  .
В этой плоскости можно задать еще одно направление через
введение единичного вектора  , а наличие двух векторов позволяет доопределить
через векторное произведение и третий единичный вектор ортонормальной системы координат  .
Так устанавливается взаимно однозначное соответствие между
расширенным до трехмерного пространством порождающего конечного элемента и
физическим или мировым декартовым пространством задачи.
Здесь  - вектор перемещения произвольной точки стержня,  -
вектор перемещения согласованной с ней точки опорной линии,  - вектор углов
поворота опорной линии, и  единичные ортогональные векторы, лежащие в плоскости
ортогональной опорной линии ,  и  проекции на эти векторы радиус-вектора от точки
опорной поверхности до средней точки сечения стержня, а  и  координаты
произвольной точки сечения стержня в базисе единичных ортонормальных векторов.
Как и для оболочек, естественно требовать, чтобы f была
кососимметричной функцией и для нее справедливы все замечания об аналогичной
функции для оболочек.
Аналогично с ситуацией для оболочки, если в некоторых точках
сечения определены перемещения, то для каждой точки сечения i справедливы
соотношения:
Как и для оболочек для удобства введем обозначение  и
запишем в матричной форме:
 => и  .
Как и в оболочках установлена линейная связь между
перемещениями всех точек сечения стержня и точки на опорной линии с
перемещениями и углами поворота точки пересечения сечения и опорной линии.
Алгоритм дальнейших действий аналогичен алгоритму, использованному при
построении матриц жесткости оболочек. Описание геометрии стержневого элемента
будем производить с помощью линейного конечного элемента, ориентированного
узлами и по программистской традиции координаты левого нижнего и правого
верхнего углов прямоугольника.
Такое представление стержня не снижает общности, так как обычно
сечение стержня описывается площадью, двумя изгибными моментами инерции и
моментом инерции на кручение, что для прямоугольника дает систему из четырех
нелинейных уравнений для четырех неизвестных, из решения которой и найдем
необходимые параметры для эквивалентного в смысле сопротивления материалов
сечения. Программа mathematica решает эту задачу, но решение оказывается слишком
громоздким, поэтому лучше использовать три параметра: площадь и два изгибных
момента инерции. Тогда решение получается достаточно простым и затем, если
полярный момент не устраивает в качестве меры кручения, это решение можно
скорректировать численным методом по четырем переменным, используя предыдущее в
качестве начального приближения. Для примера приведем решение по трем параметрам
Из множества решений естественно выбрать с положительные,
остальные решения следует опустить:
 . Для прямоугольников подкоренное выражение  тождественный
нуль.
На основании этой информации построим для каждого узла
линейного элемента базисные векторы сечения, и затем трехмерный элемент. При
этом ребра этого элемента в сечении должны иметь не менее двух узлов в каждом
направлении, а ребра в направлении образующей стержня естественно должны быть
согласованы с линейным элементом.
После этого, как и для оболочечных элементов находим матрицу
жесткости и объемные силы для элементного уравнения C U=F. Как показано ранее,
справедливы соотношения связывающие перемещения трехмерного элемента и
перемещения и углы поворота для линейного элемента U= A Q. Подставим их в
уравнения и получим C A Q=F , после этого умножим слева на транспонированную
матрицу A и придем к матрице жесткости линейного элемента и приведенным к нему
нагрузкам  .
Здесь снова по существу имеем переход к обобщенным перемещениям
и обобщенным силам, подобный обычному для аналитической механики.
Все замечания сделанные для оболочечных элементов в направлении
нормали справедливы и для линейных элементов, но уже в двух ортогональных
направлениях в сечении.
Анизотропные свойства
Компоненты тензора деформаций, используемые при записи
функционала потенциальной энергии представляются в глобальной декартовой системе
координат, поэтому при вычислениях тензор упругих констант тоже должен быть
представлен в этой же системе. В изотропном случае запись хорошо известна и не
вызывает затруднений. В любой другой ортонормальной системе координат,
полученной из глобальной системы координат смещениями и вращениями, запись будет
такой же в изотропном случае.
Для многих анизотропных сред можно выбрать такую ориентацию, в
которой тензор деформаций запишется особенно просто и обычно только в ней это
удается сделать. В такой ситуации необходимо тензор деформации, записанный в
базисе этой системы координат записать в базисе глобальной декартовой системы
координат. Подробно подобные преобразования рассматриваются в любой книжке по
тензорному исчислению. Здесь приведем лишь некоторые сведения необходимые при
построении матриц жесткости.
Сначала рассмотрим вектор X в двух системах координат -
глобальной декартовой системе с базисом (i,j,k) и локальной с единичными
направляющими векторами в качестве базиса (e1,e2,e3).
Будем скалярно умножать последовательно на i,j,k направляющие
векторы глобальной системы и объединим в вектор и получим
Вспомним, что скалярное произведение единичных векторов - это
косинус угла между ними и четная функция, поэтому направление измерения угла не
имеет значения и столбцы матрицы T по существу направляющие орты локального
базиса. На тензор в ортогональном базисе можно смотреть как на вектор,
элементами которого снова являются вектора и т.д. поэтому если меняем базис, то
надо последовательно умножать на матрицу T, поднимаясь по вложенности. В
тензорных обозначениях матрицу T можно записать  , где верхний
индекс номер строки, а нижний индекс номер столбца. Тогда если представлять
вектор в виде вложенных столбцов можно записать преобразование элементов к
глобальным координатам в тензорном виде  при этом
действует соглашение о суммировании парных диагональных индексов, то есть,
например  и все остальные. То есть здесь четыре операции суммирования. Работать
придется с большими объемами информации и монотонными действиями, поэтому
сначала подготовим несколько полезных функций.
 и 
 => .
Еще одна функция для преобразования тензора к другой ортогональной
системе координат. В матрице M столбцы это базисные векторы системы
координат, в которой описан тензор. Эти базисные векторы описаны в системе
координат, в которую надо преобразовать тензор:
После такой подготовки можно приступить к получению конкретных
соотношений. Сначала проверим хорошо известный факт независимости матрицы
упругости от системы координат в изотропном случае.
Видим, что запись матрицы упругости не изменилась. Но совсем
другая ситуация если немного поменяем исходную матрицу введением нулей:
Здесь только фрагмент реальной матрицы упругости.
Простой вид в одной системе координат сразу изменяется, и
матрица упругости становится полностью заполненной, даже при незначительных
возмущениях от изотропного закона Гука. При практических матричных вычислениях
часто учитывают симметрию тензоров напряжений и деформаций и при векторном
представлении этих тензоров используют не девятимерное, а шестимерное
пространство. При этом для сохранения возможности вычислять энергию в виде
скалярного произведения по предложению Кармана для недиагональных компонент
тензора деформаций берутся удвоенные значения компонент тензора деформаций. Но
как отмечает Работнов "нужно помнить, что представление тензора в виде вектора
имеет лишь ограниченный смысл и пригодно только для определенной фиксированной
системы отнесения". В принципе несложно получить и матрицы для преобразований
координат в шестимерном пространстве, однако, вычислительный выигрыш
незначителен, а проблем в общем случае получаем много и ясность операций уходит.
В общем случае соотношения при необходимости можно получить, используя
возможность замены координат тензора деформаций линейным оператором
сужения пространств, а потом, перенося, этот оператор на базис, в общем, вполне
аналогично тому, как осуществлялся ранее переход к обобщенным координатам на
поверхностях или линиях из трехмерной теории упругости на основе конгруэнтного
преобразования.
Практически аналогично при необходимости надо преобразовывать и
тензор упругости из косоугольной нормированной системы координат в глобальную
декартову систему координат задачи.
Есть еще несколько ситуаций, когда удобно использовать
векторное представление тензора. Ранее принимались гипотезы о поведении оболочек
и стержней для перемещений, но достаточно часто принимаются гипотезы о свойствах
тензора напряжений. Например, в оболочках и стержнях, можно принять гипотезу о
нулевых напряжениях в поперечном направлении. В тканях или сетках
достаточно разумной представляется гипотеза о нулевых сдвиговых напряжениях для
тензора напряжений, записанного в определенной системе координат. Подобных
ситуаций можно придумать достаточно много, например, при расчетах бамбуковых,
древесных, веревочных и тканевых конструкций и сооружений туристами.
Полезны эти подходы и другим экстремалам, например, при изготовлении каменных
топоров, или кремниевых микросхем.
Необходимые соотношения легко получить, преобразуя закон
Гука в векторной записи, что соответствует линейному оператору,  . Для этого
достаточно ввести в этот линейный оператор матрицу гипотез G, которая в общем
случае будет единичной, а в случае принятия гипотез о поведении компонент
тензора напряжений именно в этой матрице будет представлена соответствующая
информация  . Действуя аналогично, можем принимать и гипотезы о
поведении компонент тензора деформаций  . Логично объединить эти
выражения для того, чтобы иметь возможность принимать смешанные гипотезы  . Теперь можем записать выражение для плотности упругой энергии в виде
квадратичной формы от деформаций  , с другой стороны ничего не изменится, если в
скалярном произведении векторов поменяем сомножители местами:  . Так как
это разные формы записи одного и того же, ничего не изменится, если мы сложим и
разделим пополам один и тот же математический объект, тогда получим выражение
для плотности упругой энергии:  .
Матрица этой квадратичной формы симметричная, так как, по
существу, складываем квадратную матрицу с транспонированной матрицей. Это
представление матрицы упругих констант справедливо только в локальной системе
координат, поэтому его можно преобразовать к тензорному виду, затем найти
представление тензора в глобальной системе координат и снова перейти к
векторному представлению, более привычному для практических вычислений в методе
конечных элементов.
Приведем пример подобных преобразований на произвольных
матрицах гипотез. Построим и преобразуем к глобальной системе координат
тензор упругих констант при некоторых случайных гипотезах для
напряжений и деформаций. Получим матричное представление
этого тензора :
Получим  .
Таково представление гипотез в локальной системе координат, а в
глобальной системе:
Здесь только фрагмент матрицы упругих констант.
Видим, что упрощение представления в одной системе
координат, существенно усложняет ситуацию в другой системе координат.
Компьютер справится с этим без проблем. Попытавшись обратить эту матрицу, легко
убедимся, что она вырождена, а, следовательно, для глобальной матрицы могут
понадобиться дополнительные граничные условия, помимо обычных, вызванных тем,
что деформации, как линейные комбинации производных вносят в задачу.
Таким образом, можем рассчитывать объекты самой разной природы
и с различными свойствами. После того, как из решения задачи минимизации
квадратичного функционала найдем перемещения, можем найти и все остальные
интересующие величины. При этом иногда надо представить, например тензоры
напряжений в системе координат связанной с многообразием. При этом проще
преобразование из глобальной системы координат в локальную систему
координат осуществить в тензорном виде.
Преобразование будет происходить из локальной системы
координат в глобальную, следовательно, в формулах, подобных использованным для
преобразования тензора упругости, нам надо использовать матрицу,
транспонированную к M, так как для ортонормированных матриц транспонированная
матрица равна обратной матрице. А мы имеем дело с обратным к отображению
вектора, из локальной системы в глобальную систему координат, линейным
преобразованием.
Сначала подготовим некоторые полезные функции:
Получим  ,  ,  .
Теперь практическое преобразование в локальную систему
координат:
Получим  ,  , .
Эти преобразования можно использовать для того, чтобы
сравнивать, например, показания тензодатчиков с вычислениями, или определять
напряжения в направлении армирования, или напряжения в веревках и других
протяженных объектах.
Параметры Кармана
Объем вычислений большой, поэтому полезно
рассмотреть переход к параметрам Кармана, сокращающим объем
вычислительной работы на основе учета симметрии тензора малых
деформаций и перехода из девятимерного пространства, описывающего
деформированное состояние в малом объеме, в шестимерное пространство для
деформированного состояния.
Сначала рассмотрим логику этого преобразования. Пусть u(x) - вектор перемещения в точке
x, и пусть dx = n dx ,
где n - направление произвольного вектора.
Дифференциал вектора перемещений можем записать  , где  -
- производная Фреше или градиент вектора перемещения. Такая запись
позволяет записать ее формально в виде матричного умножения:
 => .
Обычно нас интересует изменение вдоль конкретного направления,
или другими словами - проекция на это направление, то есть скалярная
величина 
С другой стороны в силу свойств скалярного
произведения скаляр не изменится, поэтому применим к этому равенству операцию
транспонирования  .
В правой части квадратичная форма, к которой можно
применить обычную операцию симметризации, сложив эти равенства и разделив
пополам. Получим 
В силу произвольного выбора n видим, что с
помощью матрицы  можно определить относительную деформацию в любом
направлении из точки x. Записанная в виде тензора, она и назавается тензором
малых деформаций.
 =>
 .
Хотя матрица имеет 9 компонент, но в силу симмерии независимы
только 6. Следовательно, если рассматривать тензор как вектор в девятимерном
пространстве, то его можно представить как отображение из шестимерного
пространства независимых параметров.
При этом преобразовании важно проверить, чтобы
плотность упругой энергии в квадратичном функционале не изменилась. Полезно так
же понаблюдать за тем не изменяются ли при этом напряжения. А, следовательно,
можем ли говорить о сохранении закона Гука, справедливого для девятимерного
пространства, и для шестимерного. Приведенный фрагмент кода отвечает на эти
вопросы:
это показывает, что преобразование из
шестимерного в девятимерное пространство деформаций всегда корректно
осуществляется. Так выглядит матрица квадратичного функционала в нем:
 =>  - сравнение плотности упругой энергии вычисленное в девятимерном и
шестимерном пространстве совпадают.
 =>  => ,
из этих выражений следует, что для изотропных объектов и компоненты
напряжений совпадают. В изотропном случае можем говорить о сохранении закона
Гука при переходе в шестимерное пространство.
Эта благостная картина изменяется если мы
рассматриваем произвольный анизотропный случай.
матрица квадратичного функционала тоже
изменилась,  -
выражения для плотности упругой энергии тоже одинаковы в обоих размерностях, но
в касательных напряжениях видим расхождения -  и
расхождение  .
Из этого рассмотрения следует вывод о том, что
для сокращения вычислений всегда можно переходить в шестимерное пространство к
параметрам Кальмана, но при вычислении напряжений в анизотропном, по крайней
мере, случае следует использовать девятимерное пространство, или эквивалентное
тензорное представление.
Нагрузки
Поведение нагрузки можно представить в виде функции от
перемещений. При применении метода конечных элементов к теории упругости, в
качестве неизвестных параметров обычно выступают перемещения, поэтому, оставаясь
в рамках линейных задач надо ограничиться линейным приближением  , где u- вектор
перемещений,  - постоянный вектор нагрузок,  - значение
производной Фреше от перемещений по пространственным переменным в
недеформированном состоянии, по структуре квадратная матрица.
В случае перехода к обобщенным переменным для оболочек и
стержней параметры в точке будут шестимерным вектором, содержащим перемещения и
углы поворота, а матрица Фреше будет соответственно шестимерной квадратной
матрицей. В строительстве коэффициенты этой матрицы обычно называются
коэффициентами постели, потому что с их помощью обычно моделируются условия
взаимодействия сооружения с землей. Понятно, что поскольку перемещения в общем
случае в каждой точке свои, то и матрицы будут в каждой точке свои. В инженерной
практике часто неравномерностью пренебрегают и в качестве матрицы берут
осредненное значение, полученное на основе экспериментов, или фундаментальных
решений для полупространства. Такое решение для полупространства, нагруженного
прямоугольником с равномерным давлением, можно найти в книге Лурье по теории
упругости. Используя принцип суперпозиции для линейных задач, и прилагая на
нескольких разных прямоугольниках смещенных в пространстве нагрузки, всегда
можем провести достаточное количество вычислительных экспериментов для
нахождения средних значений для матриц Фреше, в дальнейшем, следуя традиции,
будем обозначать ее буквой C .
В общем случае матрица Фреше будет несимметричной, что не очень
удобно в практическом плане. Нагрузка в методе конечных элементов обычно
используется в линейном функционале, а представление ее в виде линейного
отображения от перемещений преобразует исходный функционал в сумму линейного и
квадратичного функционала. Квадратичный функционал, или другими словами
квадратичная форма порождает класс эквивалентности матриц, в качестве
представителя которого можно выбрать и симметричную матрицу, по обычному для
квадратичных форм правилу симметризации
В статье "Стационарные точки и условия на переменные" показана
возможность использования однородных или Грассмановых координат для записи
квадратичного функционала в виде  . С учетом линейной зависимости усилия от
перемещений можем записать функционал в симметричном виде  . То есть
нагрузки в этом случае модифицируют матрицу жесткости. В этой ситуации есть еще
одна особенность. Матрица A строится на основе производных, и, следовательно, с
точки зрения перемещений порождает класс эквивалентности по сдвигу на постоянный
вектор. Другими словами она строится с точностью до констант и является
полунормой, что приводит к ее вырожденности и необходимости доопределять
краевыми условиями как минимум до обратимости этой матрицы. Задание зависимости
нагрузки от перемещений делает матрицу  уже полной Соболевской нормой
в C1, или другими словами обратимой матрицей, а, следовательно, задание краевых
условий уже не становится необходимым с точки зрения обратимости и задача может
быть решена и без них. Из физических соображений, учитывая, что силы
сопротивления противоположны движению, понятно, что коэффициенты матрицы
жесткости C будут отрицательными и, следовательно,
матрица будет иметь ненулевую Соболевскую норму С1, что обеспечит возможность
разложения без дополнительных условий.
Это свойство иногда используется для автоматического
определения условий в ситуациях, когда краевые условия не важны или неизвестны,
а известна уравновешенная с точки зрения теоретической механики система сил.
Тогда матрицу A можно заменить слегка возмущенной матрицей  , с малым
параметром e при единичной
матрице, система с этой матрицей будет уже невырожденной и можно найти решение.
При этом естественно нет смысла смотреть на перемещения, а картина
деформирования будет определена с точностью до сдвигов, а напряжения будут
вполне разумными.
Как отмечает Зенкевич по поводу
относительных погрешностей в итерационных процессах, хорошая
вычислительная практика брать в качестве епсилон половину половину
машинной точности, то есть для 16 значащих цифр, будет e=10e-8. Это видимо максимальная разумная величина и
в нашем случае, реально можно брать порядка e=10e-10, и корректирровать диагональные элементы
умножением на (1+e), или более формальная запись коррекции
A+ e*diag(A).
В пользу этого вычислительного приема можно привести следующий
довод. Решение получается с некоторой точностью, то есть существует область
возмущенных матриц конечного элемента, которые дают одинаковые решения в смысле
попадания в область ограниченную точностью. Следовательно, есть область
возмущений матрицы и согласованная с ней область решений, а,
следовательно, из этой области можем выбрать такую возмущенную матрицу, которая
имеет ненулевой определитель и может быть обращена, при этом решение будет все
еще в области неразличимых между собой с точки зрения точности. Естественно,
близость может пониматься только в смысле производных, или деформаций и
напряжений. Перемещения могут быть определены с точностью до вектора
произвольного смещения как твердого тела.
Выбор возмущения логично производить на основе оценки
обусловленности системы базисных функций конечного элемента, как это
определялось в статье "Свойства изопараметрических конечных элементов".
Такие приемы полезны, например, при расчетах движения тел в
средах. Понятно, что последовательность e будет порождать последовательность
решений и уже для последовательности решений можно осуществлять предельный
переход.
Заключение
Основная цель метода конечных элементов - осуществлять
практические расчеты, а для этого нужны матрицы с хорошими свойствами. Большой
вычислительный опыт и теория обусловленности позволяют быть уверенными в
трехмерном случае. Что будет происходить с вычислительными свойствами после
описанных преобразований, столь похожих на конгруэнтные преобразования матриц?
Если возьмем на нормали две точки, то это будет конгруэнтное преобразование,
если точек больше, то можем доопределить матрицы произвольными столбцами до
квадратной матрицы, вектор с обобщенными координатами дополним нулевыми
компонентами до подходящей размерности. После этого снова получим
конгруэнтное преобразование типа
 и, следовательно, в любой ситуации оказываемся
в условиях теоремы Сильвестра, которая утверждает сохранение количества
положительных, отрицательных и нулевых собственных значений матриц в результате
конгруэнтных преобразований.
Поскольку для трехмерных элементов положительная определенность
матриц жесткости хорошо известна, и следует хотя бы из положительности
 при
положительной определенности матрицы упругости, следовательно, и матрицы
жесткости для поверхностных и линейных элементов будут положительно
определенными и не будет чрезмерного роста элементов при
 разложении.
Остается еще вопрос о том, сколько надо условий для
разрешимости матриц. Из механических соображений понятно - чтобы не было
движений как твердого тела, потому, что с точки зрения деформаций, как
комбинации частных производных, надо убрать неопределенность в матрицах
жесткости, которая возникает при дифференцировании.
Рассматривались в основном матрицы жесткости, но в методе
конечных элементов используются еще матрицы масс связанные с ускорениями и
матрицы демпфирования. Их преобразование при переходе к обобщенным координатам
осуществляется аналогично, для этого достаточно обратить внимание на тот факт,
что преобразование не зависит от времени, а зависит только от пространства.
Следовательно, соответствующие матрицы всегда можно вынести из-под знака
дифференцирования по времени или наоборот ввести под знак дифференцирования.
Существует и другой более традиционный подход к оболочкам и
стержням, связанный с использованием существующих теорий балок и оболочек. О нем
можно почитать у Зенкевича, как в классической книге, так и новых трех книгах,
изданных на английском языке. Информацию о них легко найти в Internet.
Аналогично можно поступать и в случае задач не связанных с
механикой деформируемого тела.
Литература
- Работнов Ю.Н. Механика деформируемого твердого тела: - М.:
Наука, 1988.
- Дьяконов В.П. Системы
символьной математики Mathematica 2 и Mathematica 3: - М.:СКпресс, 1998.
- Зенкевич О. Метод конечных
элементов в технике: - М.: Мир, 1975.
- Биргер И.А., Мавлютов Р.Р.
Сопротивление материалов: - М.: Наука, 1986.
- Сьярле Ф. Метод конечных
элементов для эллиптических задач: Пер. с англ. - М.:Мир, 1980.
- Васидзу К. Вариационные методы
в теории упругости и пластичности: Пер. с англ. - М: Мир,1987.
- Родионова В.А., Титаев Б.Ф.,
Черных К.Ф. Прикладная теория анизотропных пластин и оболочек: С-Пб, Изд. С-Пб. университета 1996.
- Партон В.З.,
Перлин П.И. Методы математической теории упругости:
- М.: Наука, 1981.
Санкт-Петербург, ноябрь-декабрь 2003 года
В начало |
