Full Text

Введение Рассматривается метод исследования напряженного состояния вблизи особых точек конструкций с особенностью в виде составных плоских клиньев и пространственных ребер. Поскольку вершины клиньев и ребра являются потенциальными концентраторами напряжений, изучению напряженно-деформированного состояния (НДС) в их окрестностях посвящены многочисленные публикации [1-10 и др]. Обычно авторы помещают в особую точку полюс криволинейной системы координат и сводят исследование НДС в ее окрестности к проблеме собственных чисел однородной упругой задачи. При таком подходе не учитываются реально заданные ограничения на параметры состояния непосредственно в особой точке, так как в полюсе криволинейной системы (полярной, цилиндрической или сферической) параметры состояния не определены. Более того, реально заданные ограничения в особых точках обычными методами, используемыми для исследования НДС, не могут быть учтены, потому что здесь, как правило, количество задаваемых ограничений превосходит количество задаваемых ограничений в обычной точке поверхности (контура) тела. Например, в вершине клина, образующие которого свободны от нагрузки, нулевой вектор напряжений задается на двух площадках (в обычной точке лишь на одной площадке). Из сказанного следует, что в существующих исследованиях НДС вблизи особых точек не рассматриваются ограничения, заданные непосредственно в этих точках, и имеются препятствия к учету таких ограничений, связанные с неклассическим их заданием. Данное обстоятельство обусловливает необходимость разработки новых подходов, способных строить решение в окрестности особых точек, согласующееся со всеми заданными в них ограничениями. Один из возможных подходов предлагается в настоящей работе. Он основывается на отслеживании выполнения в окрестности особой точки всех заданных в ней обязательных алгебраических равенств (ОАР), связывающих параметры состояния. В ОАР, в частности, входят граничные условия, а также условия непрерывности перемещений, напряжений и деформаций на поверхностях соединения различных элементов и т.п. Количество ОАР, формулируемых в особой точке, превышает количество ОАР в обычной точке поверхности (контура) тела. К примеру, в плоской задаче в обычной точке защемленного контура ОАР - это обращение в нуль компонентов вектора перемещений и относительного удлинения линейного элемента, направленного по касательной к контуру (всего три равенства). В точке контура, являющейся вершиной клина, образующие которого защемлены, количество ОАР равно пяти: обращение в нуль компонентов вектора перемещений относительных удлинений линейных элементов, направленных по образующим клина, и сдвига между этими элементами. Сформулированные в особой точке ОАР служат критериями достоверности решений, получаемых различными методами. Исследование ОАР как системы линейных уравнений позволяет еще до решения задачи механики [11, 12] - находить сочетания геометрических и материальных параметров конструкции, при которых особая точка теряет свой статус (перестает быть особой); - обнаруживать зависимости между геометрическими и материальными параметрами конструкции, с одной стороны, и параметрами нагрузки - с другой, обусловливающими несовместность ОАР, которая служит причиной сингулярного поведения НДС; - выявлять ограничения на параметры нагрузки, обеспечивающие корректность исследования НДС в рамках симметричной теории напряжений. В настоящей работе приводится предлагаемый авторами итерационный алгоритм и его конечно-элементная реализация. 1. Постановка задачи Рассматриваются элементы конструкции, имеющие особенности в виде плоского клина или пространственного ребра, составленные из двух различных изотропных материалов (рис. 1.1, 1.2), подвергающиеся механическому или температурному нагружению. Поверхность тела вблизи особых точек считается свободной от каких-либо внешних воздействий. Рис. 1.1 Составной клин: - внешние нормали к образующим клина; - перпендикулярные им орты Рис. 1.2. Нормальное сечение ребра Г 1.1. Плоская задача Рассматривается плоский клин (см. рис. 1.1). Составляющие элементы клина соединены по линии, касательная к которой в вершине клина принимается за ось декартовой ортонормированной системы координат . Принято: образующие клина свободны от нагрузки, углы при вершине составляющих элементов 1, 2 клина - , где (1.1) Для компонентов тензоров напряжений и деформаций приняты обозначения , , где - номер составляющего элемента; нормальные напряжения - ; касательные напряжения - . На параметры состояния в вершине А накладываются ограничения (ОАР): а) граничные условия (свободная граница) (1.2) б) условия непрерывности напряжений на линии контакта (1.3) в) условия непрерывности деформаций на линии контакта (1.4) Условия (1.4) в рассматриваемой задаче выполняются автоматически лишь в исключительных случаях. 1.2. Пространственная задача Рассматривается составное ребро (см. рис. 1.2). Под пространственным ребром Г понимается линия, образованная пересечением двух различных образующих поверхностей элемента конструкции. В точке А ребра Г построено его нормальное сечение и введен ортонормированный базис . Орт направлен во внешнюю сторону тела по касательной в точке А к линии пересечения нормальной плоскости ребра Г и поверхности соединения элементов конструкции; орт - по касательной к ребру, а так, чтобы тройка векторов была правой. С построенными ортами связана декартова система координат . Принимается, что нормали принадлежат нормальной плоскости ребра Г. Дополнительно к обозначениям п.1.1 принято: - касательные напряжения на поверхности элемента 1 ребра в направлении , - касательные напряжения на поверхности элемента 2 ребра в направлении . На параметры состояния в точках ребра Г накладываются ограничения (ОАР): а) нормальные и касательные напряжения на площадках, ориентируемых векторами , обращаются в нуль, т.е. (1.5) б) на поверхности соединения нормальные и касательные напряжения непрерывны: (1.6) в) на поверхности соединения непрерывны относительные удлинения и сдвиги: (1.7) Равенства (1.7) в данной задаче автоматически не выполняются. Эти равенства служат условиями согласования кинематических соотношений с физическими уравнениями, описывающими механическое состояние соединяемых элементов. Задача состоит в разработке алгоритма, позволяющего для подчиняющегося физическим уравнениям термоупругости составного тела построить решение, удовлетворяющее соотношениям механики деформируемого тела и ограничениям (1.2)-(1.4) в случае плоской или ограничениям (1.5)-(1.7) - в случае пространственной задачи. 2. Алгоритм построения решения 2.1. Подструктурный вариант смешанного метода конечных элементов Итерационный алгоритм построения решения на каждом шаге последовательных приближений использует метод конечных элементов (МКЭ) в форме смешанного подструктурного варианта. Применение смешанного функционала для построения разрешающих уравнений МКЭ позволяет выразить компоненты деформации в узлах КЭ-сетки (а следовательно, и напряжений) без использования операций дифференцирования приближенного решения и каких-либо методов восполнения. Функционал, в котором независимыми функциями являются перемещения и деформации, построен в работе [13], (2.1) где - матрица дифференциальных операторов для записи вектора деформаций через перемещения, т.е. ; D - матрица упругих модулей материала; - вектор температурных деформаций; W - потенциальная энергия заданных объемных и поверхностных сил. Из условий стационарности функционала (2.1) следуют уравнения равновесия, зависимости Коши и граничные условия в напряжениях. С использованием функционала (2.1) строится подструктурный вариант МКЭ. Тело V, в котором разыскивается решение, разбивается на r частей (подструктур, подобластей) так, чтобы в каждой части материальные параметры были непрерывны. Каждая часть тела (подструктура) с номером k разбивается на конечные элементы. Решение для перемещений разыскивается в классе непрерывных функций во всем теле V, а для деформаций - в классе функций, непрерывных в отдельных подструктурах. Введены параметры: - число компонент вектора перемещений в узле; - число параметров деформаций в узле; - число узлов в конечном элементе; - число узлов в подструктуре с номером k; - общее число узлов; определяют размерности векторов и матриц в конечных элементах; - размерности массивов в подструктурах, - порядок разрешающей системы линейных алгебраических уравнений МКЭ. Для рассматриваемых в работе задач (плоских и осесимметричных) выбираются четырехугольные 8-узловые конечные элементы Далее на примере плоской термоупругой задачи , формируется разрешающая система МКЭ. В конечном элементе каждой подструктуры определены: 1) вектор узловых перемещений ; 2) вектор и матрицы функций форм в локальной системе координат 3) матрица градиентов функций форм 4) вектор деформаций . В результате получаем конечно-элементные соотношения и функционал (2.1) приводится к виду (2.2) Вследствие независимости векторов перемещений и деформаций условие стационарности функционала записывается равенствами (2.3) (2.4) Для записи соотношений (2.3) и (2.4) в терминах метода конечных элементов для подструктур введены матрицы (2.5) Элементарные матрицы вычисляются с использованием формул Гаусса , , (2.6) , Здесь - квадратурные коэффициенты Гаусса в узлах интегрирования; - локальные координаты; - якобиан преобразования; - соответственно векторы узловых, поверхностных и объемных сил. Из условия (2.3) для каждой подструктуры получаем векторы деформаций (2.7) где (2.8) Если соотношение (2.7) для деформаций в подструктуре подставить в (2.4) и ввести матрицу жесткости подструктуры (2.9) получим разрешающую систему уравнений МКЭ для всей расчетной области (2.10) Размерности матриц в соотношениях (2.6)-(2.10): в конечном элементе - , в подструктуре - , , , , . Особенностью подструктурных матриц и глобальной матрицы жесткости является их сильная разреженность, что учитывается при реализации алгоритма. Решением разрешающей системы алгебраических уравнений (2.10) является глобальный вектор перемещений. Деформации в узлах подструктуры определяют соотношения (2.7). Для определения напряжений в узлах подструктуры в векторе деформаций выделяются сечения - компоненты деформаций в каждом узле, учитываются начальные деформации и производятся вычисления в соответствии с физическими уравнениями . 2.2. Системы координат, используемые в алгоритме Через X, Y обозначаются координаты точек плоского тела в глобальной ортогональной декартовой системе. В узлах КЭ-сетки, принадлежащих контуру, линии соединения элементов конструкции, и в особой точке вводятся декартовы ортогональные координаты x, y (далее узловые координаты). Ориентация узловой системы координат относительно глобальной устанавливается по следующим правилам. На рис. 2.1, а показан узел, расположенный на контуре. К сторонам смежных элементов в этом узле проведены нормали , составляющие углы и с осью X. Через обозначается угол между осями X и x. На рис. 2.1, б показан узел на линии соединения различных элементов конструкции. Как и в предыдущем случае, строятся углы и , составляемые нормалями , с осью X. Угол между осями X и x выбирается равным . На рис. 2.1, в показан узел, совпадающий с особой точкой. Ось x направляется по касательной к линии соединения составляющих элементов. Через обозначается угол между осями X и x. В случае пространственной осесимметричной задачи аналогично строятся узловые системы координат xyz. При этом плоскость xz совпадает с плоскостью RZ глобальной цилиндрической системы. а б в Рис. 2.1. Узел i: а - находится на свободной границе; б - лежит на общей линии; в - является особой точкой А Преобразование компонент напряжений при переходе от глобальных координат к узловым осуществляется с использованием матриц: - для плоской задачи - для осесимметричной (2D) задачи - для осесимметричной (3D) задачи Аналогично строится матрица преобразования компонент деформации. Эти матрицы позволяют выразить напряжения и деформации в i-м узле k-й подструктуры в узловых координатах по формулам (2.11) где - сечение, выделяемое в массиве , соответствующее этому узлу. Например, для напряжений в плоском случае будем иметь где Здесь индексы определяют компоненты - столбцы сечения , соответствующие перемещениям в i-м узле, . 2.3. ОАР в узлах КЭ-сетки Через обозначено подмножество узлов КЭ-сетки, в каждом из которых заданы какие-либо ограничения на параметры состояния. Перемещения, отвечающие узлам образуют вектор-подмножество . Ограничения в узлах в термоупругой задаче представляют собой линейные алгебраические соотношения. В зависимости от положения узла КЭ-сетки, геометрических и материальных параметров элементов конструкции алгебраические соотношения разделяются на несколько типов и подтипов. Строятся эти соотношения в узловых координатах. 2.3.1. ОАР для плоского случая 2D (ndim = 3) 1. Тип 1. В точках границы выполняются условия . 2. Тип 2. На линии соединения выполняются 3 соотношения: 3. Тип 0. Особая точка. В зависимости от геометрических и материальных параметров составляющих элементов рассмотрены следующие возможные подтипы особой точки: 3.1) коэффициенты теплового расширения на параметры состояния накладываются 6 ограничений: 3.2) параметры состояния подчиняются 6 ограничениям: 3.3) , параметры состояния подчиняются 6 ограничениям: 3.4) , на параметры состояния накладываются ограничения: Здесь 2.3.2. ОАР для осесимметричного случая 2D (ndim = 4) 1. Тип 1. В точках границы выполняются условия 2. Тип 2. На линии соединения выполняются 4 соотношения: 3. Тип 0. Особая точка. Рассмотрены следующие возможные подтипы особой точки: 3.1) и коэффициенты Пуассона , на параметры состояния накладываются 8 ограничений: 3.2) и коэффициенты Пуассона , на параметры состояния накладываются 7 ограничений: ; 3.3) , на параметры состояния накладываются 7 ограничений: 3.4) , на параметры состояния накладываются 7 ограничений: Здесь 2.4. ОАР для осесимметричного случая 3D (ndim = 6) 1. Тип 1. В точках границы выполняются условия 2. Тип 2. На линии соединения выполняются 6 соотношений 3. Тип 0. Особая точка. Рассмотрены следующие возможные подтипы особой точки. 3.1) коэффициенты Пуассона , коэффициенты теплового расширения и на параметры состояния накладываются 12 ограничений: 3.2) и - 11 ограничений: 3.3) , и - 10 ограничений: 3.4) - 11 ограничений: 3.5) , - 11 ограничений: 3.6) - 10 ограничений: Здесь Построенные ОАР с использованием равенств (2.11) в каждом случае записываются в виде матричного равенства (2.12) Матрица и вектор определяются типом (подтипом) соответствующего случая. 2.5. Процедура итерационного процесса Вектор U представляется объединением векторов меньшей размерности , а матрица - объединением двух прямоугольных матриц , в результате чего равенство (2.12) записывается в виде (2.13) Равенство (2.13) рассматривается как система уравнений для перемещений . Матрица этой системы является прямоугольной, число уравнений оказывается большим, чем число неизвестных, поэтому ее решение существует лишь в обобщенном смысле (псевдорешение). Итерационный процесс построения решения задачи организуется следующим образом: 1) перемещения считаются известными на n-1 шаге приближений, из уравнений (2.13) определяется вектор перемещений в n-м приближении; 2) перемещения рассматриваются как граничные условия в узлах при решении термоупругой задачи (2.10). В результате определяется n-е приближение для перемещений ; 3) начальное приближение находится из решения задачи, в которой в качестве ОАР учитываются лишь граничные условия. Итерационный процесс продолжается до тех пор, пока алгебраические соотношения (2.12) не будут выполняться с заданной точностью. Сходимость контролируется среднеквадратичным значением величины вектора невязок. Таким образом, процесс итерационного решения задачи проводится путем решения на каждом шаге итерации обратной задачи - поиска перемещений, обеспечивающих минимум величины вектора невязок. 3. Реализация алгоритма 3.1. Структура программного комплекса Алгоритм решения задачи реализован в виде программного комплекса, состоящего из препроцессорной подготовки данных, основной части и постпроцессорной обработки результатов. Препроцессор служит для подготовки входных данных. Входные данные для препроцессора - это конечно-элементная сетка, свойства материалов, заданные перемещения, узловые и распределенные силы, приращение температуры, которые выгружаются из программного комплекса ANSYS. Препроцессор анализирует эти данные, определяет особые точки, их тип, параметры для построения узловых координат, формирует текстовые файлы с данными о задаче и особых точках. Основная часть комплекса делится на пять блоков (вспомогательных и основных): 1) инициализация входных данных, управление процессом решения (начать решение, продолжить); 2) формирование - сборка базовой разрешающей системы МКЭ (2.10) и коэффициентов sp-подсистемы (2.13); 3) решение системы (2.13) для определения очередного приближения перемещений ; 4) решение базовой (или модифицированной на итерациях) разрешающей системы МКЭ (2.10); 5) определение параметров НДС (деформации, напряжения), оценка погрешности решения, формирование текстовых файлов для использования постпроцессором. Постпроцессор служит для визуализации результатов решения (построение изолиний, линий уровня и т.д.) [14], а также для их перезаписи в формате, пригодном для графических систем SURFER, TECPLOT. 3.2. Особенности реализации программного комплекса Описываемый программный комплекс реализован на современном языке Fortran-90/95 [15-18]. Выбор этого языка вызван тем, что на нем написано большинство высокоэффективных алгоритмов по численным методам, используемым в комплексе. В частности, алгоритм сингулярного разложения матриц для определения перемещений - решения прямоугольной sp-подсистемы [16, 17] и - разложение для решения сильноразреженной разрешающей системы (модули Ma28 из библиотеки HSL (Harwell Subroutine Library)). При исследовании полей напряжений, имеющих значительные изменения в малых областях, характерный размер конечных элементов может оказаться настолько малым, что при вычислении его площади численными методами погрешность вычислений оказывается соразмерной или даже превышающей величину этой площади. Поэтому в программном комплексе вычисления проводятся с максимально возможной в языке Fortran типом точности real_16. Другой особенностью вычислительной реализации комплекса является большое количество матричных операций. Все это обусловливает существенные затраты машинного времени при решении задач. Ускорение вычислений в программном комплексе осуществляется посредством применения технологии OpenMP и поддерживающих эту технологию компиляторов Intel, оптимально реализующих встроенные матричные и векторные функции (сложение, вычитание, скалярное произведение). Эти меры позволяют ускорить вычисления в 4-5 раз. Компиляция выполнена на суперкомпьютере с параллельной архитектурой TESLA Fermi K20 в Пермском государственном научном исследовательском университете. С использованием программного комплекса решены некоторые конкретные задачи механики деформируемого твердого тела [12, 19, 20]. 4. Напряжения при температурной нагрузке вблизи края поверхности соединения составного цилиндра (пример) В приведенном здесь примере демонстрируется сходимость предлагаемого вычислительного процесса и отличие построенного итерационного решения вблизи особой точки от решения, получаемого в классическом конечно-элементном подходе. Рис. 4.1. Составной цилиндр Рассматривается составной цилиндр высотой l = 20 мм, радиусом R= 5 мм (рис. 4.1). Элементы 1, 2, составляющие цилиндр, имеют материальные характеристики: Цилиндр подвергается однородной температурной нагрузке Обязательные алгебраические равенства на линии особых точек (граница поверхности соединения) построены в работе [19], приведены в п. 2.3.2 (подтип 3.4) и применительно к настоящему примеру, записываются десятью соотношениями: (4.1) (4.2) Ограничения (4.1) отвечают классическому подходу к решению задачи. Они представляют собой граничные условия на линии особых точек и условия непрерывности напряжений на поверхности соединения элементов 1, 2. Равенства (4.2), отражающие условия непрерывности деформаций на поверхности соединения, в классическом подходе не рассматриваются. Считается, что эти условия будут выполнены автоматически вследствие непрерывности перемещений точек элементов 1, 2 на поверхности соединения. Исследование показывает, что автоматическое выполнение условий непрерывности деформаций на поверхности соединения в особых точках возможны лишь в исключительных случаях при определенных сочетаниях материальных параметров соединяемых элементов [20]. В итерационном решении задачи использовался 4-угольный 8-узловой элемент, характерный линейный размер элемента составлял 0,1 мкм. Сходимость итерационного процесса иллюстрируется на рис. 4.2, а, где приведена зависимость среднеквадратического значения величины вектора невязки от количества итераций. В результате выполнения 380 итераций решение для напряжений, отнесенное к максимальному значению , удовлетворяет ограничениям (4.1) с погрешностью, не превышающей 0,5 %. Равенства (4.2) выполняются соответственно с точностью 0,6 и 0,3 %. Графики распределения всех напряжений в окрестности особой точки приведены в работе [19]. На рис. 4.2, б показаны напряжения (i = 1,2) вблизи особой точки, полученные на одной КЭ-сетке методом МКЭ (ANSYS) и методом итераций. Видно, что в малой окрестности точки А решения существенно различны. ANSYS-решение не удовлетворяет условию непрерывности напряжений. Условия непрерывности деформаций (4.2) в ANSYS-решении выполняются с большой погрешностью, соответственно 176 и 198 %. Вне малой окрестности особой точки решения, полученные разными методами, совпадают. 3, 4 2 1 а б Рис. 4.2. Результаты вычислений после 380 итераций: а - зависимость среднеквадратической величины вектора невязки от количества итераций; б - напряжения на линии ОА вблизи особой точки: 1, 2 - ANSYS-решение в элементах 1, 2; 3, 4 - соответствующее итерационное решение Заключение Предлагаемый в работе итерационный конечно-элементный алгоритм предназначен для построения решения упругой и термоупругой задач, согласованного со всеми заданными ограничениями непосредственно в особых точках составных плоских и пространственных конструкций. Такие решения, в частности, дают возможность оценить области достоверности результатов исследований, полученных методами механики деформируемого твердого тела, не учитывающими неклассическое задание граничных условий в особых точках. Данный подход может найти применение в исследовании полей напряжений вблизи шероховатых поверхностей, а также в изучении проблем механики композитных материалов и конструкций, в механике разрушения, в механике трещин и т.д.

About the authors

V M Pestrenin

Perm State University

I V Pestrenina

Perm State University

L V Landik

Perm State University

References

Statistics

Views

Abstract - 137

PDF (Russian) - 55