«МАРКОВ ИВАН ПЕТРОВИЧ ГРАНИЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ ОДНОРОДНЫХ ТРЕХМЕРНЫХ ЭЛЕКТРОУПРУГИХ И АНИЗОТРОПНЫХ УПРУГИХ ТЕЛ Специальность 01.02.04 – механика деформируемого твердого тела Диссертация на соискание ...»

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ОБРАЗОВАНИЯ

«НИЖЕГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМ. Н.И. ЛОБАЧЕВСКОГО»

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

На правах рукописи

МАРКОВ ИВАН ПЕТРОВИЧ

ГРАНИЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ

ОДНОРОДНЫХ ТРЕХМЕРНЫХ ЭЛЕКТРОУПРУГИХ И

АНИЗОТРОПНЫХ УПРУГИХ ТЕЛ

Специальность 01.02.04 – механика деформируемого твердого тела Диссертация на соискание ученой степени кандидата физико-математических наук

Научный руководитель:

доктор физико-математических наук, профессор Игумнов Леонид Александрович Нижний Новгород – Содержание ВВЕДЕНИЕ…………………………………………………………….…………………. Глава I. Постановки задач, метод, методика решения и программная реализация…. 1.1. Математическая модель……………………………………………………... 1.2. Гранично-элементная методика ……………..…………………………….. 1.2.1. Граничное интегральное уравнение……………………….…………. 1.2.2. Гранично-элементная дискретизация………………………………… 1.2.3. Метод квадратур сверток…………………………………………....... 1.3. Программная реализация…………………………….……………………… Глава II. Фундаментальные решения и модельные задачи равновесия……………… 2.1. Анизотропные фундаментальные решения………………………………... 2.1.1. Построение статических фундаментальных решений……………. 2.1.2. Статические функции Грина для трансверсально изотропных сред……………………………………………………………………. 2.1.3. Численные примеры построения статических функций Грина… 2.1.4. Интерполяционный подход………………………………………… 2.1.5. Динамические функции Грина, численные примеры……………. 2.2. Модельные задачи равновесия….….….….….….….….….….….….….….. 2.2.1. Однородный упругий анизотропный куб под действием нагрузки на часть торца……………..……..…..…..…..…..…..…..…..…..….. 2.2.2. Однородный электроупругий куб под действием одноосной нагрузки……………………………………………………………….. 2.2.3. Эффективность методов построения статических анизотропных фундаментальных и сингулярных решений.……………………….. Глава III. Гранично-элементное моделирование….…………………………………… 3.1. Анизотропные упругие задачи……………………….…………………….. 3.1.1. Cтатическая задача о действии давления внутри сферической полости…………………………………………………….….….…… 3.1.2. Статическая задача об анизотропном кубе с полостью………….. 3.1.3. Действие стационарной горизонтальной нагрузки на торец Г-образного однородного упругого анизотропного тела………….. 3.1.4. Одноосное стационарное растяжение упругого анизотропного призматического тела……………………………….….….….….….. 3.1.5. Действие нагрузки в виде функции Хевисайда по времени на торец Г-образного однородного упругого анизотропного тела........ 3.2. Анизотропные электроупругие задачи……………………….……………. 3.2.1. Равновесие однородного электроупругого Г-образного тела под действием разности потенциалов, приложенных к торцам.…….... действием равномерно распределенной вертикальной нагрузки 3.2.4. Контактная задача Герца для электроупругого полупространства 3.2.5. Задача об одновременном действии электрического потенциала и нагрузки в виде функции Хевисайда по времени на однородное ЗАКЛЮЧЕНИЕ…………………………………………...………………………………. СПИСОК ЛИТЕРАТУРЫ…………………………………………………...…………… Для современной техники теория упругости с сопряженными полями играет важную роль при решении возникающих задач. Многие проблемы инженерной практики сводятся к исследованию элементов конструкций, ответственных узлов и деталей машин и оборудования, а в научном плане к исследованию деформируемых тел и сред при статических или динамических нагружениях с учетом существенной анизотропии материала и связанности механических и немеханических полей.

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

С разнообразием подходов по решению динамических задач теории упругости можно познакомиться в ряде монографий [17, 21, 23, 25, 26, 59, 63, 64, 72, 73, 76, 78, 82, 102, 103, 156]. В работе развивается метод граничных элементов (МГЭ), который является распространенным универсальным численно-аналитическим методом решения широкого круга задач теории упругости с сопряженными полями. В качестве обобщающих работ по упругой статике и динамике можно привести следующие [1, 3, 4, 12-15, 20, 22, 25, 47, 61, 62, 74, 75]. Изучению процессов распространения волновых полей в средах со сложными свойствами и разработке соответствующих методов исследования посвящены следующие работы [8, 16, 22, 47, 68, 71, 96, 105, 111, 120, 123, 141]. В большинстве работ краевые и начально-краевые задачи теории упругости сводятся к интегральным уравнениям, для решения которых имеется широкий круг численных методов. С аналитическими методами решения задач динамической теории упругости можно познакомиться в монографиях [4, 18, 24, 44, 75, 83, 84].

Для решения динамических задач теории упругости большее значение имеют методы интегральных преобразований. В работе они важны в сочетании с методом граничных интегральных уравнений (ГИУ). Впервые метод интегральных преобразований при решении динамических начально-краевых задач теории упругости применен Лэмбом в 1904г. Развитие метода интегральных преобразований можно проследить по публикациям [5, 27-31, 45, 48, 67, 69, 70, 85-95, 97, 109, 120, 140, 153, 166, 177].

Относительно динамических задач анизотропной теории упругости отметим работы В.А.Свекло [79-81], В.С.Будаева [9, 10], И.Г.Филипова [101], В.А.Сарайкина и Л.И.Слепяна [77], Ю.Э.Сеницкого [86, 89, 91], А.Ф.Федечева [100] и других. Впервые МГЭ появился в работе Н.И. Мусхелишвили в 1937 г., хотя метод потенциала можно отсчитывать, например, с работ И.Ньютона.

Успех применения ГИУ обеспечен результатами, полученными в теории многомерных сингулярных интегральных уравнений в работах: С.Г. Михлина [60], В.Д.

Купрадзе, Т.Г. Гегелиа, М.О. Башелейшвили, Т.В. Бурчуладзе [98] и других. Вопросам построения ГИУ статических задач теории упругости и разработке МГЭ посвящены работы А.Я. Александрова [2], Ю.Л. Бормота (1977), Ю.В. Верюжского [11], Р.В.

Гольдштейна [19], М.И. Лазарева [49], Ю.А. Мельникова [58], О.П. Николаева [65], В.З.

Партона и П.И. Перлина [66], А.Г. Угодчикова и Н.М. Хуторянского [99], Г.И. Яха [104], T.A. Cruse [116], J.C. Lachat и J.O. Watson [50, 129], F. Paris и E. Alarcon (1980), F.J. Rizzo и D.J. Shippy [158] и многих других авторов. В ходе развития метода ГИУ и МГЭ сформировались рекомендации по использованию граничных элементов высокого порядка и для вычисления коэффициентов дискретного аналога ГИУ применять численное интегрирование с помощью формул Гаусса, а сингулярные интегралы сводить, например, к несобственным с помощью метода понижения особенности, и вычислять несобственные интегралы стандартными квадратурными формулами.

Первое решение плоских нестационарных динамических задач теории упругости было осуществлено в работах T.A. Cruse и F.J. Rizzo в 1968г. В работе G.V. Narayanan и D.E. Beskos [147] было предложено применять совместно с МГЭ для решения нестационарных динамических задач метод Дурбина [119].

Первый МГЭ-подход во временной области, был представлен в работе W.J. Mansur [136, 137]. Обобщение его подхода дано в работе H. Antes [108]. Применение МГЭ для решения произвольных двумерных нестационарных динамических задач теории упругости было осуществлено в работах J. Niwa [149] и G.D. Manolis [135]. Использование МГЭ для задач с включениями и полостями описано в работе S. Hirose [125].

Распространение традиционного МГЭ-подхода непосредственно во временной области на решение трехмерных нестационарных динамических задач теории упругости было осуществлено в работах Н.М.Хуторянского [99].

Формулировка M. Schanz и H. Antes [161], базирующаяся на методе, предложенном C. Lubich [133], позволила разрабатывать не традиционным способом применение МГЭ во временной области. Существенным преимуществом этого подхода является возможность использования для анализа волновых процессов, не зная фундаментального решения в явном времени [159, 160].

МГЭ признанно эффективен для анализа изотропных упругих конструкций.

Однако при работе с анизотропными материалами возникают известные трудности при выводе и реализации фундаментальных решений. Часто используемым методом получения фундаментальных решений в динамических задачах теории упругости является преобразование Фурье, которое в изотропном случае позволяет получить конечные выражения для фундаментальных решений в частотной области, но в анизотропном случае оно приводит к решениям, выражающимся через бесконечные интегралы. Работы V.T. Buchwald [113], M.J. Lightill [132] и M.J.P. Musgrave [146], впоследствии использованные в работе M.J.P. Musgrave [145] для сред с общей степенью анизотропии, а Payton [155] для трансверсально изотропных сред являются основополагающими в этом направлении. Также свой вклад внесли Tsvankin и Chesnokov [171], R. Burridge и ко. [114], M.N. Kazi-Aoual и др. [127] и H. Zhu [180]. Тем не менее, во всех этих работах фундаментальные решения в неявном виде требуют вычисления бесконечных интегралов и не подходят для эффективного применения в МГЭ. Интегральные представления на данный момент имеются для упругих статических анизотропных задач (J.L. Synge (1957)).

Вычисление статических фундаментальных решений с помощью интерполяционного подхода было впервые описано в работе R.D. Wilson и T.A. Cruse [179]. Полиномиальный метод построения анизотропных и электроупругих функций Грина был предложен в работах E.

Pan, F. Tonon и др. [151, 170]. Попытки разработать эффективные схемы для вычисления фундаментальных решений с использованием альтернативных аналитических форм – с использованием представления через ряды Фурье – были предприняты в работах Y.C. Shiah и др. [162-164]. C.Y. Wang и J.D. Achenbach [173, 174] подошли к задаче нахождения динамической функции Грина в анизотропных средах с помощью интегрального преобразования Радона, получив выражение для фундаментальных решений в частотной области через интегралы, определенные на поверхности единичной сферы для трехмерных задач и на единичной окружности для двухмерных задач. Статическое фундаментальное решение и представления C.Y. Wang и J.D. Achenbach дают возможность эффективной численной реализации.

В работе M. Kogl, L. Gaul (2003) [121] развит МГЭ для решения нестационарных задач теории упругости со связанными полями на основе приближенного метода с двойными использованием теоремы взаимности.

Ознакомиться с методами определения констант электроупругих материалов можно в работе И.А. Паринова и др. [106]. Первое трехмерное фундаментальное решение статического электроупругого оператора произвольной анизотропии было получено W.F.J Deeg (1980). Фундаментальные решения для динамической электроупругости были получены A.N.Norris (1994) [107, 150], а решения во временной области были получены Н.М. Хуторянским и H. Sosa (1995). Для МГЭ ценность представляет подход с представлением матрицы Грина на базе использования интегрального преобразования Радона.

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

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

Научная новизна работы заключается в следующем: гранично-элементное моделирование статических и динамических краевых/начально-краевых задач смешанного типа трехмерных линейных теорий анизотропной упругости и электроупругости реализовано на основе матриц Грина в форме Радона и шагового метода численного обращения интегрального преобразования Лапласа; в использовании комбинированного подхода к вычислению анизотропных функций Грина в точках интегрирования; в оригинальном программном обеспечении, реализующем в рамках ГЭ-подхода согласованную модель аппроксимации на основе четырехугольных элементов; в гранично-элементном решении статических анизотропных упругих и электроупругих задач; гранично-элементном моделировании стационарных задач анизотропной теории упругости; шаговом ГЭ-моделировании решения задачи о действии силы в виде функции Хевисайда по времени на торец трехмерного анизотропного упругого Г-образного однородного тела; решении нестационарной трехмерной электроупругой задачи для Г-образного тела со смешанными краевыми условиями; ГЭ-моделировании динамического прогиба композитной балки.

Достоверность исследований основана на строгом соответствии исходных начально-краевых задач в частных производных линейных теорий анизотропной упругости и электроупругости системе используемых точных граничных интегральных уравнений в изображениях по Лалпасу; на применении корректно построенных дискретных аналогов регуляризованных ГИУ; на применении оттестированных подходов к вычислению анизотропных фундаментальных решений; на тщательно проработанных алгоритмах численной реализации гранично-элементного подхода; на анализе сеточной и шаговой сходимостей получаемых ГЭ-решений и их сравнении с аналитическими результатами и с решениями других авторов.

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

Основные положения, выносимые на защиту Развитие и создание гранично-элементного методического и программного обеспечения решения краевых/начально-краевых задач статики и динамики трехмерных однородных электроупругих и анизотропных упругих тел на основе матриц Грина в форме Радона и метода квадратур сверток в форме шагового метода численного обращения преобразования Лапласа.

Гранично-элементное моделирование следующих задач равновесия:

о действии давления внутри сферической полости в анизотропном смешанная задача об упругом анизотропном кубе с полостью;

о действии разности потенциалов, приложенных к торцам однородного электроупругого Г-образного тела;

о действии равномерно распределенной вертикальной нагрузки и/или поверхностной плотности заряда на призматическое электроупругое тело;

контактная задача Герца для электроупругого полупространства;

Гранично-элементное моделирование следующих динамических задач:

о действии стационарной нагрузки на торец Г-образного однородного о действии одноосной стационарной растягивающей нагрузки на упругое о действии нагрузки в виде функции Хевисайда по времени на торец Г-образного однородного анизотропного трехмерного упругого тела;

о действии нестационарной нормальной силы на упругую композитную о действии электрического потенциала и нагрузки в виде функции Хевисайда по времени на однородное Г-образное электроупругое тело.

Результаты диссертационной работы докладывались на XII международной конференции «Современные проблемы механики сплошной среды» (Ростов-на-Дону, 2008), XIV Нижегородской сессии молодых ученых – математические науки (Н. Новгород, 2009), XV Нижегородской сессии молодых ученых – математические науки (Н. Новгород, 2010), X Всероссийском съезде по фундаментальным проблемам теоретической и прикладной механики (Н. Новгород, 2011), XVIII Международном симпозиуме «Динамические и технологические проблемы механики конструкций и сплошных сред»

им. А.Г. Горшкова (Ярополец, 2012), XIX Международном симпозиуме «Динамические и технологические проблемы механики конструкций и сплошных сред» им. А.Г. Горшкова (Ярополец, 2013), форуме молодых ученых (Н. Новгород, 2013), XXV Международной конференции «Математическое моделирование в механике деформируемых сред и конструкций. Методы граничных и конечных элементов» (Санкт-Петербург, 2013), VII Всероссийской (с международным участием) конференции по механике деформируемого твердого тела (Ростов-на-Дону, 2013), XX Международном симпозиуме «Динамические и технологические проблемы механики конструкций и сплошных сред» им. А.Г. Горшкова (Вятичи, 2014), III Международном симпозиуме «Physics and Mechanics of New Materials and Underwater Applications (PHENMA 2014)» (Кхон Кэн, Таиланд, 2014), VIII Всероссийской молодежной научно-инновационной школе (Саров, 2014).

Опубликовано 19 работ [32-43, 52-57, 126], из них 13 по теме диссертации. В изданиях, рекомендуемых ВАК РФ для опубликования результатов кандидатских диссертаций, опубликовано 6 работ в соавторстве [33, 34, 36, 42, 43, 53]. Результаты работ [33, 34, 36, 42, 43, 53] принадлежат И.П. Маркову кроме постановок задач и постпроцессорных представлений результатов исследований.

Структура и объем работы Диссертационная работа состоит из введения, трёх глав, заключения, списка литературы из 180 наименований. Общий объем диссертации составляет 168 страниц машинописного текста, включая 235 рисунков и 69 таблиц.

На различных этапах работа поддерживалась грантами Президента РФ для государственной поддержки ведущих научных школ (№ НШ-3367.2008.8 2008-2009 гг.; № НШ-4807.2010.8. 2010-2011гг.; № НШ-2843.2012.8 2012-2013гг.; № НШ-593.2014.8 2014гг.); средствами ФЦП "Научные и научно-педагогические кадры инновационной России" на 2011 - 2013 годы (ГК №14.740.11.0872 от 29 апреля 2011г., ГК №14.B37.21.1137 от 13 октября 2012г., ГК №14.B37.21.1249 от 14 сентября 2012г., ГК №14.B37.21.0471 от 3 сентября 2012г., ГК №14.B37.21.1495 от 1 октября 2012г., ГК №14.B37.21.1866 от 4 октября 2012г., ГК №14.B37.21.1902 от 4 октября 2012г.); грантами РФФИ (проект № 12-08-00984-а, № 12-08-31572-мол_а, № 12-08-00698-а, № 13-08-00531а, № 13-08-00658-а, № 14-08-00811-а, 14-08-31410-мол_а (руководитель проекта)).

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

анизотропной теории упругости, представлена математическая постановка начально-краевой задачи. Приведены уравнения движения сплошной электроупругой среды, записан закон Гаусса для случая квази-статического электрического поля.

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

Во втором параграфе записан вид обобщенной начально-краевой задачи после применения интегрального преобразования Лапласа по временной переменной при отсутствии объемных сил и нулевых начальных условиях. На основе теоремы взаимности Бетти и формулы Сомильяны приведены ГИУ прямого подхода в пространстве изображений по Лапласу. Представлена гранично-элементная дискретизация на основе регуляризованных ГИУ. Приведена гранично-элементная методика решения ГИУ на основе согласованной поэлементной аппроксимации, позволяющая получить решения, параметризованные комплексным параметром преобразования Лапласа. Кратко описан метод квадратур сверток для численного обращения интегрального преобразования Лапласа. Его частный случай в форме шагового метода используется для получения решения во временной области.

гранично-элементного моделирования трехмерных краевых задач анизотропной упругости и электроупругости на базе матриц Грина в форме Радона и метода квадратур сверток.

фундаментальных решений (функций Грина) для линейных теорий упругости и электроупругости. Описаны различные подходы к вычислению статической части функций Грина. Приведены выражения статических функций Грина для трансверсально изотропных упругих и электроупругих сред. Проведено тестирование численной реализации подходов к построению статических фундаментальных решений. Приведены полутоновые трехмерные визуализации статических функций Грина на единичной сфере для упругих материалов с различной степенью анизотропии и для полностью анизотропного электроупругого материала. Проведен ряд численных экспериментов для верификации численной реализации схемы вычисления динамических анизотропных фундаментальных решений.

фундаментальных решений в пространстве изображений по Лапласу в виде суммы сингулярной и регулярной частей для упругих и электроупругих сред, полученные на основе метода с применением интегрального преобразования Радона. Описаны следующие подходы к построению статических анизотропных упругих и электроупругих функций Грина: интегральный, полиномиальный, с применением рядов Фурье и интерполяционный. Приведенные выражения статических функций Грина для трансверсально изотропных упругих и электроупругих сред использовались для тестирования численной реализации методов построения статических фундаментальных решений. Представлены трехмерные визуализации в виде поверхностей компонент анизотропной упругой статической матрицы Грина в точках на единичной сфере для четырех различных анизотропных материалов: монокристалла алюминия, ниобата бария-натрия, сапфира и углепластика. Для электроупругой среды визуализации приведены для полностью анизотропного материала. Для установления работоспособности и корректности численной реализации схемы построения динамических анизотропных фундаментальных решений проведен ряд численных экспериментов.

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

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

трехмерных задач статики и динамики анизотропных и электроупругих однородных тел.

В первом параграфе решены следующие статические упругие задачи: о действии давления внутри сферической полости, расположенной в неограниченной анизотропной упругой среде; смешанная задача о полости внутри анизотропного упругого куба. Рассмотрены две стационарные задачи: о действии нагрузки на торец однородного анизотропного Г-образного упругого тела; об одноосном растяжении упругого анизотропного призматического тела. Полученные на разных ГЭ-сетках результаты позволяют сделать вывод о наличии сеточной сходимости; ГЭ-решения приведены в сравнении с аналитическими решениями и результатами других авторов. Решена задача анизотропной упругой динамики о действии нагрузки в виде функции Хевисайда по времени на торец Г-образного однородного анизотропного упругого тела. Даны сравнения полученных ГЭ-решений с конечно-элементными (КЭ) и ГЭ-решениями других авторов, полученных по приближенному методу с двойным использованием теоремы взаимности. Помимо сеточной сходимости исследована сходимость по параметру шагового метода обращения преобразования Лапласа. Представлены решения задачи о динамическом изгибе композитной балки от действия нестационарной силы. ГЭ-решения сравнивались с данными по динамическому эксперименту и КЭ-решениями, полученными в программных комплексах конечно-элементного моделирования «Динамика-3»* и ANSYS. Продемонстрировано хорошее соответствие полученных численных решений результатам эксперимента.

Второй параграф посвящен гранично-элементному моделированию задач электроупругости. Решены следующие статические задачи: о действии разности потенциалов, приложенной к Г-образному однородному электроупругому телу; о действии вертикальной силы или заданной поверхностной плотности заряда на торец призматического электроупругого тела; о действии нагрузки на дневную поверхность электроупругого полупространства; контактная задача Герца для электроупругого полупространства и контактная задача Герца с дефектом. Проведено сравнение полученных гранично-элементных решений с известными аналитическими решениями и результатами других авторов. Дано решение задачи электроупругой динамики о совместном действии нагрузки в виде функции Хевисайда по времени и электрического потенциала на Г-образное однородное электроупругое тело. Исследована сеточная и шаговая сходимость. Приведено сравнение ГЭ-решений с КЭ- и ГЭ-решениями других авторов, полученных по приближенному методу с двойным применением теоремы взаимности.

диссертационном исследовании.

Сертификат соответствия Госстандарта России № POCC RU.ME20.H Глава I Постановки задач, метод, методика решения и программная реализация 1.1. Математическая модель Основными соотношениями линейной анизотропной теории упругости являются уравнение движения сплошной среды (уравнение Коши), выражающее баланс импульса:

и закон сохранения импульса:

Эти уравнения дополняются выражениями для линейного тензора малых деформаций Коши-Грина и обобщенным законом Гука:

где ij – компоненты тензора напряжений, f j – компоненты вектора объемной силы, – плотность материала, ui – компоненты вектора перемещений, Cijkl – тензор модулей упругости и ij – компоненты линейного тензора деформации.

Компоненты вектора поверхностных усилий определяются как:

Объединяя эти уравнения, можно получить уравнения движения (уравнения Навье) для анизотропного упругого тела R 3 с границей Г:

Чтобы полностью поставить начально-краевую задачу, уравнения движения необходимо дополнить начальными и граничными условиями:

Рассмотрим однородную линейную электроупругую среду. При допущении, что электрическое поле является квази-статическим, уравнения движения сплошной электроупругой среды и закон Гаусса имеют следующий вид [115, 121, 154, 168]:

где ij – компоненты тензора напряжений, f j – компоненты вектора объемной силы, – плотность материала, ui – компоненты вектора упругих перемещений, Di – компоненты вектора электрической индукции, – плотность электрического заряда. В силу того, что пьезоэлектрические материалы являются диэлектриками, они не содержат свободных электрических зарядов: 0.

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

где Cijkl – тензор модулей упругости, измеренный при постоянном электрическом поле;

напряженности электрического поля; eijk – компоненты пьезоэлектрического тензора, ij – компоненты тензора диэлектрической проницаемости, измеренные при постоянной деформации.

Квази-статическое электрическое поле является потенциальным:

где – электрический потенциал. Тензор модулей упругости Cijkl, пьезоэлектрический тензор eijk и тензор диэлектрической проницаемости ij удовлетворяют следующим условиям симметрии:

таким образом, для случая общей анизотропии, они содержат по 21, 18 и 6 независимых материальных констант, соответственно. Кроме того, тензоры Cijkl и ij являются положительно определенными:

для любых ненулевых вещественных aij и bi.

С использованием определения линейного тензора деформации связанные дифференциальные уравнения движения для линейной электроупругой принимают вид:

Компоненты вектора поверхностных усилий и поверхностная плотность электрического заряда определяются как:

где n j – компоненты вектора единичной нормали к площадке.

Введем сокращенные обозначения для объединения электрических и упругих переменных в обобщенные векторы и тензоры [121]. Тензоры обобщенных напряжений и деформации примут вид:

поверхностных усилий и обобщенной объемной силы:

электроупругости:

Из симметрии материальных тензоров следует, что обобщенный тензор Cijkl так же удовлетворяет условиям симметрии:

таким образом, в случае общей анизотропии электроупругого материала тензор Cijkl содержит до 45 независимых материальных констант.

Обобщенный закон Гука и выражение для вектора обобщенных поверхностных усилий примут вид:

Начально-краевая задача для электроупругого тела R 3 с границей Г запишется в виде:

где jk – дельта Кронекера.

Граничные условия:

Начальные условия:

начально-краевую задачу для анизотропного упругого или электроупругого тела R 3 с границей Г:

Граничные условия:

Начальные условия:

где N – размерность задачи; N 3 для упругости и N 4 в случае электроупругости.

1.2. Гранично-элементная методика динамической упругости и электроупругости решаются методом граничных элементов.

1.2.1. Граничное интегральное уравнение В качестве метода решения описанной начально-краевой задачи в обобщенном виде используется метод граничных интегральных уравнений (ГИУ) в пространстве изображений по Лапласу.

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

где u k и t k – компоненты векторов обобщенных перемещений и усилий, соответственно, s i, 0 – параметр преобразования Лапласа. В дальнейшем предполагается, что все выкладки приводятся в пространстве изображений по Лапласу.

Пусть состояние твердого тела R 3 описывается парой переменных состояния заполненное такой же средой, что и область, и характеризуемое переменными состояния ( ij, ij ). В силу симметрии обобщенного тензора Cijkl справедлива теорема взаимности:

Соответственно, можно записать и формулу Сомильяны [99]:

где g ij x,y, s – матрица фундаментальных решений (функции Грина) уравнения движения, h jm (x, y, s) Сijkl g km,l ni, где ni – компоненты вектора единичной нормали к поверхности Г в точке интегрирования x, y – точка коллокации.

Используя формулу Сомильяны возможно определить значения обобщенных перемещений только во внутренних точках области для заданного значения параметра преобразования Лапласа s. Чтобы найти неизвестные граничные функции необходимо построить граничное интегральное уравнение (ГИУ).

Если в формуле Сомильяны осуществить предельный переход точки y из области на границу Г, то получим следующее ГИУ:

в смысле главного значения по Коши.

Полученное ГИУ справедливо как для упругих, так и для электроупругих сред.

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

1.2.2. Гранично-элементная дискретизация Для получения дискретного представления ГИУ используется согласованный регулярного представления ГИУ:

(сингулярная) часть матрицы Неймана h jk.

элементов:

где M – число граничных элементов.

В трехмерном МГЭ граничные элементы являются элементами двумерной поверхности. Для вычисления интегралов возникающих в отображается на некий контрольный элемент с локальной системой координат (1, 2 ).

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

где N e – число узлов на элементе; (e, n) – глобальный номер узла, имеющего в e-ом элементе локальный номер n.

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

Координаты 1 и 2 упорядочены так, что направление вектора n единичной нормали к границе Г является внешним по отношению к области (рис. 1).

В рамках согласованной модели аппроксимация обобщенных перемещений u осуществляется билинейными элементами, а обобщенных поверхностных сил t постоянными элементами. Для фиксированного значения параметра интегрального преобразования Лапласа s имеют место следующие выражения для обобщенных граничных перемещений и обобщенных поверхностных сил внутри элемента e :

где Rn (1, 2 ) – линейные функции формы для линейного четырехугольного элемента:

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

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

Применение метода коллокации на узлах y m аппроксимации граничных функций позволяет получить дискретный аналог ГИУ. В итоге сформируются системы линейных алгебраических уравнений (СЛАУ):

Уравнения записаны соответственно в узлах аппроксимации обобщенных граничных перемещений и обобщенных поверхностных сил.

После соответствующих преобразований получаем разрешающую СЛАУ:

где [A] – полностью заполненная несимметричная матрица, {X} – вектор неизвестных величин, в котором объединены все соответствующие компоненты искомых обобщенных граничных функций, {B} – вектор правой части.

Задача из двух подобластей (кусочно-однородная) решается добавлением еще одного ГИУ со всем комплексом гранично-элементных аппроксимаций и сведением ГИУ к системе линейных алгебраических уравнений. Системы линейных алгебраических уравнений, записанные для каждой из подобластей, при использовании условий жесткого контакта позволяют составить итоговую систему линейных алгебраических уравнений.

Коэффициентами разрешающей СЛАУ являются интегралы по граничным элементам от ядер ГИУ с некоторыми весовыми функциями, которые определяются аппроксимацией граничных функций. В зависимости от взаимного расположения точки коллокации и граничного элемента, по которому производится интегрирование, можно выделить два типа интегралов:

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

сингулярный интеграл: точка коллокации находится на элементе. В этом случае подынтегральное выражение становится сингулярным для зависимости от типа фундаментального решения различают слабосингулярные (особенность вида 1 r ) и сильносингулярные (особенность вида 1 r ) интегралы.

Для вычисления слабосингулярных интегралов применяется преобразование ЛашаВатсона, при котором элемент подразбивается на два треугольных (вырожденных четырехугольных) элемента. В случае сильносингулярных интегралов применяется метод многоуровневых иерархических подразбиений [5].

уравнений параметризованы комплексным параметром интегрального преобразования Лапласа s. Для получения решения в явном времени требуется подходящая схема численного обращения преобразования Лапласа.

1.2.3. Метод квадратур сверток следующими соотношениями:

где s i – комплексная переменная, расположенная в полуплоскости 0.

Рассмотрим метод, применяемый для численного обращения преобразования Лапласа, опирающийся на теорему о свертках, – метод квадратур сверток.

Заменим интеграл свертки (1.1) квадратурной суммой, весовые множители которой определяются при помощи изображения по Лапласу f и линейного многошагового метода [133, 134, 161]:

где n (t ) представимы численно; на основе формулы трапеций с постоянным шагом 2 / L получим:

Аппроксимация при выводе формул (1.2)-(1.3), основана на использовании линейного многошагового метода с характеристической функцией (z ) для решения задачи Коши дифференциального уравнения первого порядка, возникающей в процессе преобразования (1.1):

Многошаговый метод должен иметь порядок точности p 1, являясь строго нульустойчивым или A-устойчивым. Функция f (s) должна быть ограничена в правой полуплоскости относительно прямой c i, c i :

При условии аналитичности f (s) и ограниченности в области где, критерий устойчивости можно ослабить до A( ) устойчивости.

В качестве соответствующих примеров многошаговых методов можно привести методы дифференцирования назад порядка p 6 : для A-устойчивого метода второго порядка ( 90 ) можно записать:

Пусть функция f (s) в уравнении (1.3) вычисляется с погрешностью, тогда выбор L N и R n допускает погрешность вычисления n порядка ( ).

Метод квадратур сверток близок к шаговому методу обращения интегрального преобразования Лапласа. Но если метод квадратур сверток основан на теореме о свертке двух оригиналов, то шаговый метод обращения преобразования Лапласа основан на теореме об интегрировании оригинала. Шаговый метод можно рассматривать как частный случай метода квадратур сверток при g ( ) 1. В этом случае имеем:

В работе применяется шаговый метод обращения преобразования Лапласа при условии f (t ) y (t ).

1.3. Программная реализация Для численной реализации представленной гранично-элементной схемы был разработан программный комплекс для расчета статических и динамических задач трехмерных анизотропных теорий упругости и электроупругости. Комплекс состоит из трех независимых модулей: вспомогательной программы ANISGN, расчетного модуля ABEM и программы STINV. Все программы написаны на алгоритмическом языке Фортран и представляют собой консольные приложения. Программа ANISGN является препроцессором и генерирует интерполяционные таблицы для заданного анизотропного материала. Расчетный модуль ABEM реализует численное решение поставленной задачи в изображениях по Лапласу. Программа STINV используется для обращения полученного решения во временную область.

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

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

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

SYST GMAT GEON BCON FREQ INTG

В модуле SYST производится инициализация банка данных, затем в него считывается конфигурационный файл. Считывание параметров материала из вспомогательного файла и дальнейшая их обработка производится в модуле GMAT.

Построение полной дискретной модели геометрии задачи происходит в модуле GEON. В модуле BCON производится обработка граничных условий задачи. Модуль FREQ считывает набор значений параметра преобразования Лапласа из вспомогательного файла.

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

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

Для успешной эксплуатации программного комплекса необходимо подготовить конфигурационный файл и полный набор вспомогательных файлов. Структура конфигурационного файла организована по логической структуре программы ABEM.

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

В качестве основных выводов по главе I можно сформулировать следующее.

Глава I содержит постановки начально-краевых задач линейных теорий анизотропной упругости и электроупругости. С помощью введенных сокращенных обозначений – векторов обобщенных перемещений и усилий, сформулирована обобщенная начально-краевая задача, позволяющая рассматривать задачи анизотропной теории упругости и электроупругости с единых позиций. В пространстве изображений по Лапласу записан обобщенный вид начально-краевой задачи. Дана схема сведения начально-краевой задачи в обобщенном виде к ГИУ прямого подхода на основе теоремы взаимности Бетти и формулы Сомильяны. Представлена гранично-элементная дискретизация на основе регуляризованных ГИУ. Приведена гранично-элементная методика решения ГИУ в изображениях по Лапласу на основе ГЭ-схемы с согласованной поэлементной аппроксимацией. Для получения решения во временной области используется частный случай метода квадратур сверток – шаговый метод численного обращения интегрального преобразования Лапласа. Дано краткое описание созданного программного комплекса компьютерного моделирования трехмерных начально-краевых задач анизотропной упругости и электроупругости с использованием метода ГИУ и интегрального преобразования Лапласа. Комплекс программ состоит из трех независимых компонент. В главе описано назначение каждой из этих компонент. Указан принцип взаимодействия и структура обмена данными между этими компонентами на базе вспомогательных файлов.

Глава II Фундаментальные решения и модельные задачи равновесия 2.1. Анизотропные фундаментальные решения Фундаментальные решения или функции Грина играют ключевую роль в методе граничных интегральных уравнений и методе граничных элементов для анализа статических и динамических задач линейной теории упругости с сопряженными полями.

Рассмотрим неограниченную однородную анизотропную линейную упругую или электроупругую среду. Обобщенные динамические функции Грина Gkp (x, t ) (матрица N N, где N 3 для упругости и N 4 для электроупругости) во временной области математически определяются следующей системой дифференциальных уравнений в частных производных:

с начальными условиями:

где jp – дельта Кронекера, (x) – дельта-функция Дирака, – плотность среды, x ( x1, x2, x3 ) – точка наблюдения, t – время.

Для материалов с наиболее общей степенью анизотропии, функции Грина невозможно получить в явном конечном виде. C помощью интегрального преобразования Радона [110, 122], для фундаментальных решений и их производных возможно получить интегральные представления [121]. Пусть f ( x ) R 3, p – действительное число, тогда прямое и обратное преобразования Радона имеет вид:

Обобщенные динамические функции Грина в пространстве изображений по Лапласу с параметром s имеют вид [150, 174]:

g S – статическая (сингулярная), g R – динамическая (регулярная) части обобщенного фундаментального решения.

В упругом случае динамическая часть имеет вид [138, 150, 174]:

Для электроупругости динамическая часть имеет вид [176]:

где где m, E pm в упругом случае – собственные числа и соответствующие собственные вектора матрицы Кристоффеля jk (n) Cijkl ni nl ; в электроупругом случае – собственные числа и соответствующие собственные вектора матрицы L jp, Q – количество различных собственных чисел; cm, k m – фазовые скорости и волновые числа; n – единичный вектор направления распространения волны; D R – область интегрирования динамической части решения (половина единичной сферы).

Статическая часть матрицы фундаментальных решений (или статическая матрица Грина) для упругости и электроупругости имеет вид [121, 174]:

где D S – область интегрирования статической части (единичная окружность) 2.1.1. Построение статических фундаментальных решений В выражении (2.2) для статической части обобщенного анизотропного фундаментального решения можно заметить, что матрица Кристоффеля является симметричной и положительно определенной, поэтому обратная к ней матрица всегда существует. Учитывая это, последнее выражение в (2.2) можно напрямую вычислять при помощи подходящей квадратуры. Такой подход к вычислению статической части фундаментального решения называют интегральным [121, 174].

Статические функции Грина, определенные в (2.2), можно записать как:

где A jk (n) – присоединенная матрица для jk (n), D(n) - определитель jk (n).

Тогда, с использованием основной теоремы о вычетах, фундаментальные решения можно выразить в виде [151, 170]:

где v – произвольный единичный вектор, отличный от e, Im m 0, m 1, N – корни многочлена степени 2N: D(p q) 0, коэффициент при старшей степени. Единственным численным шагом при вычислении статической части функции Грина по формуле (2.3) является нахождение корней многочлена. Данный подход называется полиномиальным [151, 170].

Если перейти к сферической системе координат (r,, ), то статическую часть матрицы Грина можно выразить следующим образом [164, 169]:

где тензор H S (, ) является периодическим с интервалом 2 для обоих сферических углов. Соответственно, его можно представить в виде разложения в двойной ряд Фурье:

Коэффициенты Фурье (2.6) могут быть вычислены при помощи, например, квадратуры Гаусса-Лежандра. Как только коэффициенты Фурье определены, статические обобщенные функции Грина в произвольной точке можно записать как:

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

Такой метод вычисления фундаментальных решений называют методом с использованием рядов Фурье [164, 169].

2.1.2. Статические функции Грина для трансверсально изотропных сред Не считая хорошо известного решения Кельвина-Томсона для изотропных материалов, явные выражения для статических функций Грина были получены только для трансверсально изотропных сред Lifshitz and Rozenweigх [131]; Krner [128]; Willis [178];

Pan and Chou [152]. Различными авторами предпринимались попытки получить статические фундаментальные решения в явном виде для материалов с кубической симметрией, но удалось вывести только приближенные формулировки Mura and Kinoshita [144]. Дальнейшие исследования ориентированы на использование численных методов Mura [143] Tonon et al. [170].

Рассмотрим важный случай – трансверсально изотропную среду. Широко используемые материалы, например, пьезокерамические композиты и некоторые поликристаллы, относятся именно к этому классу материалов. Для трансверсально изотропных материалов тензор модулей упругости имеет следующий вид:

Приведем выражения для статических функций Грина в виде, в котором они представлены в работе Ahmad Pouya [157]. Не теряя общности, в дальнейшем предполагается, что источник расположен в начале координат, а x( x1, x2, x3 ) – точка наблюдения. Введем следующие обозначения:

Если c c13 2c44 0 :

Если c c13 2c44 0 :

Дополнительно введем вспомогательные функции:

Статические функции Грина для упругой трансверсально изотропной среды запишутся в виде:

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

электроупругих сред были получены в работе [139]. Приведем эти выражения. Для этого введем следующие обозначения:

где si – корни кубического уравнения:

a c11 (11c33 2e15e33 ) 11c13 (c13 2c44 ) c44 ( 33c11 e31 ) 2e15c13 (e31 e15 ), Функции Грина для трансверсально изотропного электроупругого материала запишутся в виде:

где Выражения для A2 j ), A3 j ) и D2, D3 получаются с помощью циклической перестановки индексов 1, 2 и 3 в выражениях для A1( j ) и D1, соответственно.

2.1.3. Численные примеры построения статических функций Грина Чтобы установить работоспособность и точность численной реализации описанных выше подходов к вычислению обобщенных анизотропных статических функций Грина рассмотрим несколько примеров. Пусть задана линейная упругая трансверсально изотропная среда, в которой плоскость изотропии параллельна координатной плоскости x1 x2. Тензор модулей упругости материала имеет следующий вид [170]:

фундаментальными решениями для трансверсально изотропной среды представлены в таблице 1, а их производные – в таблице 2.

Для интегрального метода в таблицах 3, 4 приведены результаты при различных порядках k квадратуры Гаусса-Лежандра.

1,1 4.01415886104E-07 4.01415414110E-07 4.01415886177E-07 4.01415886104E- 1,2; 2,1 -2.93155291433E-08 -2.92921654507E-08 -2.93155291563E-08 -2.93155291433E- 1,3; 3,1 -2.15178851716E-07 -2.15191821967E-07 -2.15178852003E-07 -2.15178851716E- 2 2 3.88223897990E-07 3.88233939657E-07 3.88223898057E-07 3.88223897990E- 2,3; 3,2 1.72143081372E-07 1.72153457573E-07 1.72143081602E-07 1.72143081372E- 3 3 1.90032202843E-06 1.90026882712E-06 1.90032202654E-06 1.90032202843E- интегрального подхода уже при использовании простой квадратуры Гаусса-Лежандра с узлами.

Чтобы удостовериться в корректности численной реализации метода основанного на рядах Фурье, рассмотрим полностью анизотропный материал, со следующим тензором модулей упругости [164]:

В таблицах 5, 6 приведены результаты при различных значениях (параметр метода) для точки наблюдения x (r 1.0, 4, 3) в сравнении с результатами из [164].

Рассмотрим вычисление статической матрицы Грина и её производных для электроупругой среды. В качестве материала взят трансверсально изотропный пьезокерамик цирконат-титанат свинца PZT-4, со следующими параметрами [151]:

Точка наблюдения расположена в x (1, 1, 1). Результаты вычислений для компонент статической матрицы Грина и их производных по полиномиальному методу приведены в таблицах 7, 8, по интегральному методу в таблицах 9, 10, по методу с использованием рядов Фурье в таблицах 11, 12, в сравнении с результатами, полученными по явным формулам для электроупругих трансверсально изотропных материалов.

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

2.1.4. Интерполяционный подход Опишем подход к построению обобщенных анизотропных статических фундаментальных решений на основе интерполяционной схемы вычисления. Если вернуться к представлению функций Грина в сферической системе координат (2.4):

можно заметить, что тензор H S (, ) полностью определен на области изменения сферических углов и. В силу того, что область определения сферических углов ограничена, значение g S в произвольной точке может быть получено при помощи подходящей интерполяционной схемы [121, 179] для вычисления значений тензора H S, например, линейной интерполяции Лагранжа.

В качестве примера приведем трехмерные визуализации в виде поверхностей (рис. 3-26) для компонент анизотропной упругой статической матрицы Грина в точках наблюдения, расположенных на единичной сфере (r 1,, 0 ), для четырех различных анизотропных материалов со следующими тензорами модулей упругости [121, 175]:

Рис. 3 g 11 для монокристалла алюминия Рис. 4 g 22 для монокристалла алюминия Рис. 5 g 33 для монокристалла алюминия Рис. 6 g 12 для монокристалла алюминия Рис. 7 g 13 для монокристалла алюминия Рис. 9 g 11 для ниобата бария-натрия Рис. 10 g 22 для ниобата бария-натрия Рис. 11 g 33 для ниобата бария-натрия Рис. 12 g 12 для ниобата бария-натрия Рис. 13 g 13 для ниобата бария-натрия Рис. 14 g 23 для ниобата бария-натрия Рис. 19 g 13 для сапфира Рис. 20 g 23 для сапфира Рис. 21 g 11 для углепластика Рис. 23 g 33 для углепластика В случае электроупругой среды рассмотрим полностью анизотропный материал со следующими параметрами:

Рис. Рис. Рис. Рис. Рис. Рис. Рис. Рис. Рис. 2.1.5. Динамические функции Грина, численные примеры Для достижения высокой точности вычисления динамических функций Грина и Неймана применяется составная квадратурная формула Гаусса-Лежандра [172]. Пусть n – количество участков разбиения промежутка интегрирования [a, b], p – количество узлов на каждом из участков, тогда:

где B j, t j – абсциссы и узлы обычной квадратурной формулы Гаусса порядка p.

Для оценки достоверности и аккуратности численной реализации схемы вычисления анизотропного динамического фундаментального решения по формуле (2.1) были выполнены три теста [138]:

Тест 1 включает в себя вычисление изотропных динамических матриц Грина и Неймана, для которых известны выражения в явном виде.

в тесте 2 производится сравнение матриц Грина и Неймана для трансверсально изотропного в плоскости x1 x2 материала, с соответствующими результатами из Тест 3 содержит сравнение матриц Грина и Неймана для полностью анизотропного материла, полученного из трансверсального изотропного путем преобразования системы координат. А именно, в декартовой системе координат {xi }; i 1, компоненты тензора модулей упругости для трансверсально изотропного в плоскости x1 x2 материала определяются следующим образом:

Поворачиваем исходную систему координат {xi } в новую {xi}, соответствующим преобразованием:

где lij1) – компоненты матрицы преобразования поворота вокруг произвольного единичного вектора v на некоторый угол.

В новой системе координат {xi} тензор модулей упругости будет иметь вид:

Дополнительное вращение на угол вокруг некоторого вектора v, заданного уже в новой системе координат:

порождает полностью анизотропный материал со следующим тензором модулей упругости:

или Функции Грина для полностью анизотропного материала g ik могут быть выражены через функции Грина для трансверсально изотропного материала g pk следующим образом:

Статическая часть матрицы Грина во всех тестах вычислялась по интегральному методу с порядком квадратуры Гаусса k 48.

Тест 1. В случае, когда материал является изотропным, численная схема должна давать результаты, близкие к получаемым по известным аналитическим выражениям.

Рассмотрим вычисление функций Грина и Неймана для случая со следующими параметрами:

где r – радиус-вектор точки наблюдения (точка нагружения расположена в начале координат), s – значение параметра преобразования Лапласа, модуль упругости и модуль сдвига модельного изотропного материала, n(n1, n2, n3 ) – единичный вектор, необходимый для вычисления матрицы Неймана hij Cmjkl gik,l nm, n – количество участков разбиения промежутка интегрирования, p – количество узлов на каждом из участков.

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

здесь и далее g обозначает матрицу Грина, полученную по аналитической формуле (тест 1 и тест 2) или с помощью соответствующего преобразования системы координат (тест 3), В таблицах 13, 14 приведено сравнение полученных результатов для функций Грина и Неймана, соответственно. Погрешность вычислений для матрицы Грина в сравнении с аналитическим результатом составила err(g, g) 5.35 10 -15, для матрицы Неймана err(h, h) 1.26 10 -15.

1,1 2.76401753E-02 -3.61239597E-02 2.76401753E-02 -3.61239597E- 1,2 1.25675660E-02 -2.84706901E-03 1.25675660E-02 -2.84706901E- 1,3 1.02613747E-02 -2.32462212E-03 1.02613747E-02 -2.32462212E- 2,1 1.25675660E-02 -2.84706901E-03 1.25675660E-02 -2.84706901E- 2,2 2.76401753E-02 -3.61239597E-02 2.76401753E-02 -3.61239597E- 2,3 1.02613747E-02 -2.32462212E-03 1.02613747E-02 -2.32462212E- 3,1 1.02613747E-02 -2.32462212E-03 1.02613747E-02 -2.32462212E- 3,2 1.02613747E-02 -2.32462212E-03 1.02613747E-02 -2.32462212E- 3,3 2.34509867E-02 -3.51749367E-02 2.34509867E-02 -3.51749367E- 1,1 -9.69358646E-02 2.49588819E-02 -9.69358646E-02 2.49588819E- 1,2 -7.30495146E-02 7.58020596E-03 -7.30495146E-02 7.58020596E- 1,3 -5.80710265E-02 5.45414213E-03 -5.80710265E-02 5.45414213E- 2,1 -7.30495146E-02 7.58020596E-03 -7.30495146E-02 7.58020596E- 2,2 -9.69358646E-02 2.49588819E-02 -9.69358646E-02 2.49588819E- 2,3 -5.80710265E-02 5.45414213E-03 -5.80710265E-02 5.45414213E- 3,1 -6.12009481E-02 7.32148647E-03 -6.12009481E-02 7.32148647E- 3,2 -6.12009481E-02 7.32148647E-03 -6.12009481E-02 7.32148647E- 3,3 -7.25718331E-02 2.27564624E-02 -7.25718331E-02 2.27564624E- Результаты теста 1 демонстрируют высокую степень точности используемой схемы численного интегрирования при вычислении как полей перемещений, так и полей поверхностных усилий.

Тест 2. Рассматривается вычисление матрицы Грина для трансверсально изотропного материала, для следующего набора параметров:

Значения компонент матрицы Грина в сравнении с результатами Dravinski и Zheng [118] приведены в таблице 15. Тест показывает высокую степень согласованности полученных результатов с результатами других авторов. Также вычисления проводились на широком наборе точек наблюдения 1 r a 15. Здесь и далее a 1 м – параметр расстояния. Изображения реальных и мнимых частей некоторых компонент матриц Грина и Неймана приведены на рис. 37-42.

1,1 3.305547E-02 7.758555E-02 3.30554788E-02 7.75855491E- 1,2 2.975494E-02 8.06068E-02 2.97549312E-02 8.06067600E- 1,3 1.88524E-03 1.9895E-04 1.88524919E-03 1.98951656E- 2,1 2.975494E-02 8.06068E-03 2.97549312E-02 8.06067600E- 2,2 3.305547E-02 7.758555E-02 3.30554788E-02 7.75855491E- 2,3 1.88524E-03 1.9895E-04 1.88524919E-03 1.98951656E- 3,1 1.88524E-03 1.9895E-04 1.88524919E-03 1.98951656E- 3,2 1.88524E-03 1.9895E-04 1.88524919E-03 1.98951656E- 3,3 1.248244E-02 1.540302E-02 1.24819052E-02 1.54030226E- g g g Тест 3. В последнем тесте сравниваются матрица Грина для полностью анизотропного материала, полученная через преобразование системы координат с матрицей, прямо вычисленной по формуле (2.1). Параметры преобразования поворота системы координат выбраны следующие:

Тензор модулей упругости трансверсально изотропного материала из теста 2 в новой системе координат примет вид [138]:

В случае следующих параметров вычисления:

результаты теста для матрицы Грина приведены в таблице 16.

1,1 4.17386762E-02 7.84042759E-02 4.17386768E-02 7.84042759E- 1,2 2.74444349E-02 7.75930501E-03 2.74444346E-02 7.75930501E- 1,3 6.42909601E-03 1.09334814E-02 6.42909356E-03 1.09334814E- 2,1 2.74444349E-02 7.75930501E-03 2.74444346E-02 7.75930501E- 2,2 2.31071190E-02 7.48237442E-02 2.31071189E-02 7.48237442E- 2,3 5.99822328E-03 -3.07601451E-04 5.99822513E-03 -3.07601451E- 3,1 6.42909601E-03 1.09334814E-02 6.42909356E-03 1.09334814E- 3,2 5.99822328E-03 -3.07601451E-04 5.99822513E-03 -3.07601451E- 3,3 1.37470675E-02 1.73461007E-02 1.37470736E-02 1.73461007E- Результаты вычислений для набора точек наблюдения 1 r / a 15 приведены на рис. 43, 44. Для полностью анизотропного материала численная реализация схемы вычисления матрицы Грина по формуле (2.1) позволяет получать результаты с высокой степенью точности.

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

2.2. Модельные задачи равновесия 2.2.1. Однородный упругий анизотропный куб под действием нагрузки на часть торца [42, 121] фундаментальных решений на точность вычислений, рассмотрим статическую упругую анизотропную задачу, конфигурация которой изображена на рис. 45. Нижний торец x3 0 м куба жестко закреплен, на 1/5 части правой x1 0.01 м стороны приложена равномерно распределенная нагрузка ~ t 0, t0 1108 Па, направленная вдоль оси x1.

Остальная поверхность тела является свободной от поверхностных усилий. Анизотропные фундаментальные решения в точках интегрирования вычислялись с помощью линейной интерполяции Лагранжа на сетках 100100, 200 200, 400 400, и 800 800. В узлах интерполяции фундаментальные решения вычислялись по интегральному методу с порядком квадратуры Гаусса равным 48. Кроме того, в случае изотропного материала фундаментальные решения вычислялись по формулам Кельвина.

изотропная сталь, трансверсально-изотропный цирконат-титанат свинца (ЦТС), ортотропная древесина и моноклинный углепластик. Тензоры модулей упругости материалов приведены в таблице 17 [121].

Расчеты проводились на трех ГЭ-сетках, сведения о которых представлены в таблице 18. Сетка «а» изображена на рис. 46, сетка «б» на рис. 47, сетка «в» - на рис. 48.

полученных решений, при использовании различных интерполяционных сеток для вычисления фундаментальны решений, с ГЭ-решениями из [121] приведено в таблицах 19-22 для стали, ЦТС, древесины и углепластика, соответственно.

Из полученных результатов можно сделать вывод, что с увеличением степени анизотропии материала, для аккуратной интерполяции фундаментальных решений требуется увеличение размера интерполяционной сетки. Как правило, сетка позволяет получить достаточно точные результаты для перемещений.

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

Единичный электроупругий куб, с центром в начале координат, подвергнут одноосному растяжению t3 100 Па вдоль ось x3, или действию равномерно распределённой поверхностной плотности заряда q 1010 Кл м 2 (рис 49). В узлах, расположенных на серединной линии x3 0.5 м задан электрический потенциал 0 В.

Остальная поверхность куба является свободной от поверхностных усилий и поверхностной плотности заряда. Анизотропные фундаментальные решения в точках интегрирования строились с помощью линейной интерполяции Лагранжа на сетках трансверсально-изотропных материалов. В узлах интерполяции фундаментальные решения вычислялись по интегральному методу с порядком квадратуры Гаусса равным 48.

пьезокерамика – цирконат-титанат свинца PZT-4, со следующими параметрами [124]:

и цирконат-титанат свинца PZT-5H, со следующими параметрами [34]:

Ось трансверсальной изотропии материалов совпадает с осью x3 используемой декартовой системы координат. В силу симметрии задачи гранично-элементная сетка строится с учетом трех плоскостей симметрии. Расчеты проводились с использованием ГЭ-сеток с различной степенью дискретизации на 1/8 части куба. Сведения о них представлены в таблице 23.

Маркер ГЭ-сетки Количество элементов Количество узлов Рассматриваются отклики перемещений u1, u3 и электрического потенциала в точке 0.5, 0.5, 0.5м. Для материала PZT-4 в таблицах 24-26 дано сравнение полученных при механической нагрузке гранично-элементных решений, в зависимости от уровня дискретизации сетки интерполяции фундаментальных решений. При нагрузке в виде поверхностной плотности заряда сравнения представлены в таблицах 27-29.

Сетка интерполяции Аналитические Сетка интерполяции Аналитические формулы Сетка интерполяции Аналитические формулы Сетка интерполяции Аналитические формулы Сетка интерполяции Аналитические формулы Сетка интерполяции Аналитические Полученные результаты свидетельствуют о сходимости гранично-элементного решения как по количеству элементов в ГЭ-сетке, так и по уровню дискретизации интерполяционной сетки фундаментальных решений. Для материала PZT-4 в таблице дано сравнение ГЭ решений при механической нагрузке, полученных на сетке интерполяции 800 800, с ГЭ- и аналитическими решениями из [124], для нагрузки в виде поверхностной плотности заряда – в таблице 31. Соответствующие сравнения для PZT-5H представлены в таблицах 32, 33.

ГЭ решение на сетке «в» -0.11631264 0.41186913 1. Аналитическое решение из [124] ГЭ решение на сетке «в» -0.59028666 0.13212334 -0. Аналитическое решение из [124] Аналитическое решение из [124] Аналитическое решение из [124] Эффективность методов построения статических анизотропных фундаментальных и сингулярных решений.

Чтобы сравнить вычислительную эффективность представленных методов построения статических анизотропных фундаментальных и сингулярных решений, при непосредственном применении в разработанном гранично-элементном программном комплексе, рассмотрим модельную статическую задачу упругости. Куб со стороной 1 м.

подвергнут действию равномерно распределенной растягивающей нагрузки ~ t 0, t0 1011 Па, приложенной на торце x3 1 м и действующей в направлении оси x (рис. 50). На противоположном торце задана скользящая заделка x3 0 м. В качестве анизотропного материала взят альфа-кварц со следующим тензором модулей упругости [167]:

Задача решалась на четырех ГЭ-сетках с различной степенью дискретизации.

Сведения о них представлены в таблице 34.

Для методов построения фундаментальных решений были выбраны следующие параметры. В интегральном подходе использовалась квадратура Гаусса с числом узлов равным 48. В подходе через ряды Фурье параметр метода был выбран равным 20. В интерполяционном подходе использовалась сетка 800 800. В таблице 35 дано сравнение аналитического решения из [51] с ГЭ-решениями для перемещений u 2 в точке A 0.5, 0.5, 0.5 м, полученными всеми четырьмя методами на всех ГЭ-сетках.

Интегральный -0.068125910 -0.068573379 -0.068732356 -0. Полиномиальный -0.068125861 -0.068573783 -0.068732893 -0. Интерполяционный -0.068125443 -0.068567061 -0.068725716 -0. Аналитическое решение Все четыре метода демонстрируют сеточную сходимость к аналитическому решению и дают очень близкие результаты. В таблице 36 приведены длины временных промежутков (в секундах), которые потребовались для полного решения задачи при использовании каждого из четырех методов.

Полиномиальный 13.907842 122.904028 1125.197593 4448. Совокупность полученных результатов позволяет сделать однозначный вывод, что с практической точки зрения интерполяционный подход к построению фундаментальных и сингулярных анизотропных статических решений, имеет наибольшую вычислительную эффективность. В случае серьезных ограничений по объему доступной оперативной памяти имеет смысл использовать интегральный подход.

В качестве итогов второй главы можно привести следующие результаты.

Представлены выражения для анизотропных динамических фундаментальных решений для упругих и электроупругих сред в пространстве изображений по Лапласу, записанные в виде суммы статической и динамической частей, полученные с помощью применения интегрального преобразования Радона. Описаны четыре метода вычисления статических анизотропных упругих и электроупругих функций Грина: интегральный, полиномиальный, с применением рядов Фурье и интерполяционный. Даны выражения статических функций Грина для трансверсально изотропных упругих и электроупругих сред. Проведено тестирование численной реализации подходов к построению статических фундаментальных решений. Представлены трехмерные полутоновые визуализации в виде поверхностей для компонент анизотропной упругой статической матрицы Грина для четырех различных анизотропных материалов. В случае электроупругой среды визуализации даны для полностью анизотропного материала. Проведена верификация численной реализации схемы вычисления динамических анизотропных фундаментальных решений. Проведен ряд численных экспериментов. Исследована работоспособность и вычислительная эффективность описанных методов построения статических анизотропных фундаментальных и сингулярных решений, при непосредственном применении в разработанном гранично-элементном программном комплексе. Решены две модельные задачи для верификации предложенной гранично-элементной методики решения статических задач линейных теорий анизотропной упругости и электроупругости с интерполяционным способом вычисления функций Грина. Полученные ГЭ-решения приведены в сравнении с аналитическими результатами и результатами других авторов.

Глава III Гранично-элементное моделирование 3.1. Анизотропные упругие задачи 3.1.1. Статическая задача о действии давления внутри сферической Рассмотрим сферическую полость с радиусом R 1 м, заданную в неограниченной анизотропной упругой среде (рис. 51). К границе полости приложено равномерно распределенное нормальное давление P 107 Па. Исследуются перемещения на границе полости. Анизотропные фундаментальные решения в точках интегрирования строились при помощи интерполяционного подхода на сетке 800 800. В узлах интерполяции фундаментальные решения вычислялись интегральным методом с порядком квадратуры Гаусса равным 48.

Решение задачи рассматривалось для четырех различных материалов: изотропный поликристаллический алюминий, монокристалл алюминия, ниобат бария-натрия и сапфир. Тензоры модулей упругости материалов приведены в таблице 37 [175].

Поликристаллический алюминий Монокристалл алюминия Расчеты проводились на трех гранично-элементных сетках с различной степенью дискретизации. Сведения о них представлены в таблице 38. Сетка «а» изображена на рис.

52, сетка «б» на рис. 53, сетка «в» - на рис. 54.

расположенных на границе полости в плоскостях x1, x2, 0 и 0, x2, x3. Для изотропного материала сравнение полученных гранично-элементных решений в точках на плоскости x1 x2 с аналитическим и ГЭ-решением из [175] приведено на рис. 55, где arctan x2 x1.

Для монокристалла алюминия сравнения в точках на плоскостях x1 x2 и x2 x3 даны соответственно на рис. 56 и рис. 57, где arctanx2 x3, для ниобата бария-натрия – на рис. 58, 59 и для сапфира – на рис. 60, 61.

ur, м Рис. 55 Перемещения u r в плоскости x1 x2 для поликристаллического алюминия ur, м Рис. 56 Перемещения u r в плоскости x1 x2 для монокристалла алюминия ur, м Рис. 57 Перемещения u r в плоскости x2 x3 для монокристалла алюминия ur, м ur, м ur, м ur, м Продемонстрировано поведение радиальных перемещений в двух координатных плоскостях на границе сферической полости для четырех материалов с различной степенью анизотропии. Проводя сравнение полученных ГЭ-решений с аналитическим решением в изотропном случае, необходимо отметить, что результаты на ГЭ-сетке «в»

значительно лучше, чем результаты на ГЭ-сетке «а», что свидетельствует о сходимости к аналитическому решению. Сеточная сходимость гранично-элементного решения наблюдается для всех рассмотренных анизотропных материалов.

3.1.2. Статическая задача об анизотропном кубе с полостью [167] Рассмотрим сферическую полость с радиусом r расположенную внутри анизотропного упругого куба с длиной ребра L 40 м, при этом стороны куба ориентированы параллельно координатным плоскостям декартовой системы координат.

На трех сторонах куба задана скользящая заделка, а на противоположных – нормальная поверхностные усилия величиной t 0 1 Па (рис. 62).

Рассматривался полностью анизотропный материал со следующим тензором модулей упругости [167]:

Расчеты проводились на гранично-элементной сетке, содержащей 600 элементов – 384 элемента на сторонах куба и 216 на полости. Исследовались абсолютные перемещения ua u12 u2 u3 в узлах, расположенных на границе полости в плоскости x1, x2, 0. Задача рассматривалась для четырех значений величины радиуса полости конечно-элементного анализа ANSYS (лицензия ANSYS Academic Research, Customer Number 623640) было получено конечно-элементное решение задачи. КЭ-сетка содержит 74612 элементов типа SOLID186 и 102564 узла. Сравнения полученных КЭ- и ГЭрешений для нормализованных абсолютных перемещений c11ua R 0 t 0 приведены на рис.

63-66, где arctanx2 x1.

Рис. Рис. Рис. Приведены абсолютные перемещения в узлах, расположенных на границе сферической полости, находящейся внутри анизотропного упругого куба. Проведено сравнение полученных ГЭ- и КЭ-решений для различных значения величины радиуса полости.

3.1.3. Действие стационарной горизонтальной нагрузки на торец Г-образного однородного упругого анизотропного тела [121] анизотропного тела, изображенного на рис. 67. На стороне x1 0.1 м приложены поверхностные усилия ~ t 0, t0 1108 Па, направленные вдоль оси x1. На нижнем торце тело жестко закреплено, остальная поверхность является свободной от поверхностных усилий. Задача рассматривалась для двух материалов: изотропной стали (с параметрами материала: E 1.008 1011 Па, 0.3 и плотностью 7800 кг м3 ), и трансверсально-изотропного цирконата-титаната свинца (ЦТС) PZT с плотностью 7850 кг м3 и следующим тензором модулей упругости [121]:

Параметр преобразования Лапласа, для которого проводились расчеты, в случае стали равен s 0 i 100, в случае ЦТС s 0 i 50. Задача решалась на двух ГЭ-сетках, сведения о которых представлены в таблице 39. Сетка «а» изображена на рис. 68, сетка «б» – на рис. 69.

Исследовался отклик перемещений u1 в узлах вдоль линии Сравнение полученных ГЭ-решений с решениями из [121] приведено для стали на рис. и для ЦТС на рис. 71. На рис. 72-73 представлены трехмерные полутоновые визуализации распределения поля перемещений u1 для стали и ЦТС, соответственно.

u1, м u1, м Приведены отклики перемещений в узлах вдоль вертикальной линии на границе тела. Полученные ГЭ-решения для стали и ЦТС на двух ГЭ-сетках даны в сравнении с результатами других авторов. Представлены трехмерные полутоновые визуализации распределения полей перемещений на сетке «б».

Рассмотрим упругую анизотропную стационарную задачу, изображенную на рис. 74. Нижний торец жестко закреплен, к верхнему торцу приложены поверхностные усилия ~ t, t 100 Па, направленные вдоль оси x.

В качестве материала взят сапфир с плотностью 3970 кг м3 и следующим тензором модулей упругости [142]:

Расчеты проводились для двух различных значений параметра преобразования Лапласа: s1 0 i и s2 0 i 4. Задача решалась на двух ГЭ-сетках, сведения о которых представлены в таблице 40. Сетка «а» изображена на рис. 75, сетка «б» – на рис.

76.

0.2, 0.15, x3 м. Сравнения полученных ГЭ-решений с решениями из [142] даны на рис.

77-79 для s1 и на рис. 80-82 для s2. На рис. 83, 84 представлены трехмерные полутоновые визуализации распределения поля абсолютных перемещений ua u12 u2 u3 для s1 и s2, соответственно.

u1, м Рис. 77 Перемещения u1 при параметре преобразования Лапласа s u2, м Рис. 78 Перемещения u 2 при параметре преобразования Лапласа s u3, м Рис. 79 Перемещения u 3 при параметре преобразования Лапласа s u1, м Рис. 80 Перемещения u1 при параметре преобразования Лапласа s u2, м u3, м Рис. 83 Перемещения u a для параметра преобразования Лапласа s преобразования Лапласа. Дано сравнение полученных на различных ГЭ-сетках решений с результатами других авторов. Представлены полутоновые визуализации распределения поля абсолютных перемещений для обоих значений параметра преобразования Лапласа.

3.1.5. Действие нагрузки в виде функции Хевисайда по времени на торец Г-образного однородного упругого анизотропного тела [36, 121] Рассмотрим задачу о действии равномерно распределённой горизонтальной нагрузки в виде функции Хевисайда по времени ~ (t ) t0 H (t ), t0 1105 Па на боковую поверхность однородного упругого анизотропного тела, изображенного на рис. 85.

Расчеты проводились на четырех ГЭ-сетках, сведения о которых представлены в таблице 41. Сетка «а» изображена на рис. 86, сетка «б» – на рис. 87, сетка «в» – на рис. 88, сетка «г» – на рис. 89.

Задача рассматривается для двух материалов с различной степенью анизотропии:

трансверсально изотропный цирконат-титанат свинца (материал 1) и моноклинный углепластик (материал 2). Плотности (в кг/м3) и ненулевые компоненты тензоров модулей упругости (в ГПа) для этих материалов приведены в таблице 42 [121]. Для идентификации графиков откликов перемещений и поверхностных усилий, построенных на сетке «а»

используется маркер «», на сетке «б» – «», на сетке «в» – «», на сетке «г» – «*».

Рассмотрим точку A с координатами (0, 0, 0.1)м. Сходимость решения в перемещениях для материала 1 продемонстрирована на рис. 90-92, для материала 2 на рис. 93-95.

Значение параметра шагового метода численного обращения интегрального преобразования Лапласа выбрано t 5 10 6 с.

u1, м u2, м u3, м u1, м u2, м u3, м Рассмотрим центральный элемент на закрепленном торце. Сеточная сходимость решения для поверхностных усилий для материала 1 продемонстрирована на рис. 96-98, для материала 2 на рис. 99-101. Рассматривались только три сетки – «а», «б» и «г», так как на сетке «в» центральный элемент отсутствует.

t1, Па t2, Па t3, Па t1, Па t2, Па t3, Па преобразования Лапласа по параметру t, продемонстрирована для материала 2 на рис.

102-104 на примере перемещений в т. А (0, 0, 0.1)м. В качестве расчетной ГЭ-сетки здесь и далее используется сетка «в».

u1, м u2, м u3, м На рис. 105, 106 представлено сравнение полученных ГЭ-решений для перемещений по оси x1 в точке А (0, 0, 0.1)м с решениями из [121] для метериалов 1 и 2, соответственно.

u1, м u1, м На рис. 107-118 приведена полутоновая визуализация поля абсолютных перемещений решения задачи для материала 1, причем палитра моделируется линейкой от 0 м (синий цвет) до 6.3 10 6 м (красный цвет). На рис. 119-130 приведена трехмерная визуализация решения задачи для материала 2. В этом случае, для наглядности, Представленные изображения соответствуют интервалу времени t [0,0.0005] с с шагом t 0,000045 c.

Рис. 109 t 0.00009 с Рис. 110 t 0.000135 с Рис. 111 t 0.00018 с Рис. 112 t 0.000225 с Рис. 113 t 0.00027 с Рис. 114 t 0.000315 с Рис. 115 t 0.000360 с Рис. 116 t 0.000405 с Рис. 117 t 0.00045 с Рис. 118 t 0.0005 с Рис. 121 t 0.00009 с Рис. 122 t 0.000135 с Рис. 123 t 0.00018 с Рис. 124 t 0.000225 с Рис. 125 t 0.00027 с Рис. 126 t 0.000315 с Рис. 127 t 0.000360 с Рис. 128 t 0.000405 с Рис. 129 t 0.00045 с Рис. 130 t 0.0005 с Для оценки влияния степени анизотропии материала на отклик перемещений, на рис. 131, 132 представлено сравнение ГЭ-решений u1 и u3 в точке B (0.05, 0.025, 0.05)м для обоих материалов.

u1, м u3, м Приведены динамические отклики перемещений для двух материалов с различной степенью анизотропии. Исследована сеточная сходимость и сходимость по параметру шагового метода численного обращения интегрального преобразования Лапласа. Дано сравнение полученных ГЭ-решений с ГЭ- и КЭ-решениями других авторов. Приведены полутоновые визуализации распространения поля абсолютных перемещений для материала 1 и трехмерные визуализации решения задачи для материала 2.

3.1.6 Динамический изгиб композитной балки Рассмотрим задачу о численном моделировании динамического испытания на изгиб композитной балки (рис. 133). Размеры образца составляют 6 мм 15 мм 80 мм. На одной стороне на полосе шириной a 1.4 мм, расположенной симметрично относительно t0 1.074748571108 Па, вид функции f (t ) дан на рис. 134. На противоположной стороне относительно начала координат. Расстояние между опорами составляет c 45 мм. На этих опорах задана скользящая заделка u1 0. Остальная поверхность балки свободна от поверхностных усилий. Использовался ортотропный материал.

Гранично-элементные решения сравнивались с результатами эксперимента (методика описана в [6, 7, 46, 112]) и соответствующими результатами, полученными в ANSYS (лицензия ANSYS Academic Research, Customer Number 623640). Расчеты проводились на двух гранично-элементных сетках с различной степенью дискретизации.

Сведения о них представлены в таблице 43. ГЭ-сетка «а» изображена на рис. 135, ГЭсетка «б» – на рис. 136. Значение параметра шагового метода численного обращения интегрального преобразования Лапласа выбрано t 1 10 7 с. Исследовался отклик перемещений u1 (t ) в точке O(0, 0, 0).

Сертификат соответствия Госстандарта России № POCC RU.ME20.H осуществлялось в системе «нагружающий стержень – опорная трубка». Схема представлена на рис. 137.

Мерные трубка и стержень были изготовлены из алюминиевого сплава Д16.

Физические и геометрические характеристики экспериментальной установки приведены в таблице 44. Фотография экспериментальной установки дана на рис. 138.

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

где I, R T – падающий, отраженный и прошедший импульсы деформации, соответственно.

Сравнение полученных ГЭ-решений прогиба балки с экспериментальным данными дано на рис. 139. Сравнения с КЭ-решениями из ANSYS и «Динамика-3»

представлены на рис. 140, 141, соответственно.

u1, м u1, м u1, м Приведено гранично-элементное моделирование динамического испытания на изгиб композитной балки. Продемонстрирована хорошая согласованность полученных на двух различных сетках ГЭ-решений и результатов эксперимента. Кроме того дано сравнение результатов с решениями, полученными в программных комплексах конечно-элементного моделирования.

3.2. Анизотропные электроупругие задачи 3.2.1. Равновесие однородного электроупругого Г-образного тела под действием разности потенциалов, приложенных к торцам [34, 121] Рассматривается статическая электроупругая задача о действии разности потенциалов, приложенных к однородному телу (рис. 142). На нижнем торце x3 0 м тело жестко закреплено и к нему приложен электрический потенциал 1 0 В. На правом торце x1 0.1 м приложен электрический потенциал 2 100 кВ. Остальная поверхность тела является свободной от поверхностных усилий и поверхностной плотности заряда.

Задача рассматривалась для трансверсально изотропного материала – цирконата-титаната свинца (ЦТС) PZT со следующими параметрами [121]:

Расчеты проводились на двух ГЭ-сетках, сведения о которых представлены в таблице 45. Сетка «а» изображена на рис. 143, сетка «б» – на рис. 144.

Исследовался отклик перемещений u3 и электрического потенциала в точках вдоль линии 0.05, 0, x3 м. Сравнение полученных решений с ГЭ-решениями из [121] приведено на рис. 145 для u3 и на рис. 146 для. На рис. 147, 148 представлены электрического потенциала, соответственно.

u3, м Приведены отклики перемещений и электрического потенциала. Проведено сравнение полученных ГЭ-решений с ГЭ-решениями других авторов. Представлены электрического потенциала.

3.2.2. Равновесие призматического электроупругого тела, под действием равномерно распределенной вертикальной нагрузки и/или поверхностной Рассмотрим гранично-элементное решение статической задачи электроупругости, когда призматическое тело подвергнуто действию равномерно распределенной механической, электрической или смешанной нагрузок. У единичного электроупругого куба на нижнем торце x3 0.5 м задана скользящая заделка – u3 0 м, 0 В. В первом варианте нагружения к верхнему торцу x3 0.5 м приложена равномерно распределённая растягивающая нагрузка ~ P 1108 Па, q 0 Кл м 2 (рис. 149).

Остальная поверхность тела является свободной – на ней заданы поверхностные усилия ti 0 Па (i 1,3) и поверхностная плотность заряда q 0 Кл м.

В качестве материала рассматриваются два трансверсально изотропных пьезокерамика – цирконат-титанат свинца PZT-4, со следующими параметрами [117]:

и цирконат-титанат свинца PZT-5H, со следующими параметрами [117]:

Ось трансверсальной изотропии материалов совпадает с осью x3 используемой декартовой системы координат. В силу симметрии задачи гранично-элементная сетка строится с учетом двух плоскостей симметрии. Расчеты проводились с использованием ГЭ-сеток с различной степенью дискретизации. Сведения о них представлены в таблице 46. Сетка «а» изображена на рис. 150, сетка «б» на рис. 151, сетка «в» на рис. 152.

Маркер ГЭ-сетки Количество элементов Количество узлов Рассматривались гранично-элементные решения задачи в точках A (0.5, 0, 0)м, B (-0.25, 0.5, 0.25)м и O (0, 0, 0)м. Статические фундаментальные решения строились двумя способами: с помощью интерполяции на сетке 800 800, значения функций Грина и её производных в узлах интерполяции вычислялись по интегральному методу с порядком квадратуры Гаусса равным 48; по аналитическим формулам для трансверсально изотропных материалов. Для задачи существует аналитическое решение [148], которое используется для сравнения с полученными ГЭ-решениями. В таблицах 47-49 приведено сравнение полученных гранично-элементных решений задачи для материала PZT-4 с численными и аналитическими результатами из [148] для точек A, B и O, соответственно.

Для PZT-5H результаты приведены в таблицах 50-52.

Точка A Точка B Точка O Точка A Точка B Точка O перемещений и электрического потенциала в точке B, в зависимости от количества граничных элементов на ГЭ-сетке. Сравнение проводилось с аналитическими решениями из [148]. Максимальная относительная погрешность составила не многим более 0.5% на ГЭ-сетке с наименьшим числом элементов.

Относительная погрешность, % ГЭ-решений с численными и аналитическими результатами из [148] приведены в таблицах 53-55 для PZT-4 и в таблицах 56-58 для PZT-5H.

Точка A Точка B Точка O Точка A Точка B Точка O в точке B, в зависимости от количества граничных элементов на ГЭ-сетке отображены на рис. 155, 156. Сравнение проводилось с аналитическими решениями из [148].

Максимальная относительная погрешность в случае электрической нагрузки получилась на порядок меньше чем в случае механической, и составила чуть более 0.045%.

Относительная погрешность, % Относительная погрешность, % Рассмотрим случай смешанной нагрузки. На верхнем торце тела одновременно заданы растягивающая нагрузка t3 P 1 108 Па и поверхностная плотность заряда аналитическими результатами из [148] приведены в таблицах 59-61 для PZT-4 и в таблицах 62-64 для PZT-5H.

Точка A Точка B Точка O Точка A Точка B Точка O являются суперпозицией значений, полученных отдельно при механической и электрической нагрузках. Зависимости относительных погрешностей от количества граничных элементов изображены на рис. 157, 158 для PZT-4 и PZT-5H, соответственно.

Максимальная относительная погрешность в случае смешанной нагрузки составила 0.323%.

Относительная погрешность, % перемещений ua u12 u2 u3 и электрического потенциала соответственно.

Механическая Электрическая Смешанная Механическая Электрическая Смешанная Рассмотрим решение задачи, когда материал PZT-4 повернут относительно рабочей системы координат на 18 вокруг оси x3, затем на 35.25 вокруг оси x1 и на вокруг оси x2. В таком случае материал становится полностью анизотропным со следующими параметрами [148]:

Нижний торец тела жестко закреплен и на нем задан нулевой электрический потенциал.

Боковые поверхности свободны от усилий и поверхностной плотности заряда. На верхнем торце приложена нагрузка ~ P 1108 Па, q 0 Кл м 2. На рис. 159 дана полутоновая визуализация распределения поля абсолютных перемещений, на рис. 160 – электрического потенциала. В силу полной анизотропии используемого материала, перемещения и электрический потенциал в направлении действии нагрузки не являются линейными.

-0. 3.2.3. Задача о действии нагрузки на дневную поверхность электроупругого полупространства [130].

Рассмотрим гранично-элементные решения статической начально-краевой задачи [130] о действии равномерно распределённой вертикальной нагрузки ~ P 121 10 9 Па на часть дневной поверхности электроупругого полупространства (рис. 161). Дневная поверхность является свободной – на ней заданы поверхностные усилия ti 0 Па (i 1,3) и ограниченного окружностью радиусом r a 1 м с центром в начале координат.

В качестве материала взят трансверсально изотропный цирконат-титанат свинца – пьезокерамик PZT5A со следующими параметрами [130]:

Исследовались отклики некоторых компонент вектора обобщенных перемещений и вектора напряжений, вычисленные по формуле Сомильяны, во внутренних точках полупространства A (1, 0, -2)м, B (2, 0, -2)м, C (3, 0, -2)м, E (4, 0, -2)м, расположенных на глубине x3 2 м вдоль оси x1. Гранично-элементная сетка строится с учетом одной плоскости симметрии. Половина расчетной сетки содержит 696 элементов и 741 узел (рис.

162).

В таблице 67 приведены полученные ГЭ-решения для нормализованных поверхностных усилий t3 P, в таблице 68 – поверхностной плотности заряда qc11 e33P в сравнении с ГЭ-, КЭ- и аналитическими решениями из [130].

перемещений u1, u3 и электрического потенциала в узлах вдоль оси x1 на дневной поверхности и величин откликов u1, u3 и, вычисленных по формуле Сомильяны, для узлов расположенных вдоль оси x1 на глубине d 2 м от дневной поверхности.

Отметим, что нормализованные перемещения u3 и электрический потенциал имеют идентичную форму отклика.

Рассмотрим вариант задачи когда исследуемые точки расположены на глубине d 1 м от дневной поверхности. Гранично-элементная сетка строится с учетом двух плоскостей симметрии. Четверть расчетной сетки содержит 348 элементов и узлов (рис. 166). На рис. 167-170 продемонстрировано сравнение полученных по формуле Сомильяны решений для нормализованных перемещений u3c11 aP, электрического потенциала 10e33 aP, поверхностной плотности заряда qc11 e33 P и напряжений t3 P соответствующими ГЭ- и КЭ-решениями из [130].

u3c11/aP 10e33/aP -qc11/e33P -t3/P Рассмотрим вариант задачи, когда исследуемые точки расположены на оси x3, а в качестве нагрузки на части дневной поверхности задана поверхностная плотность заряда q Q 15.8 Кл м 2 (рис. 171). Гранично-элементная сетка строится с учетом двух плоскостей симметрии. Четверть расчетной сетки содержит 348 элементов и 387 узлов (рис. 172). На рис. 173-175 продемонстрировано сравнение полученных по формуле координатами из [130].

-10tre33/c11Q u3e33/aQ -e33/c11aQ Решена задача о действии нормальной силы на часть дневной поверхности электроупругого полупространства. Рассмотрено два различных случая расположения исследуемых точек внутри полупространства. Решен вариант задачи, когда на части дневной поверхности задана поверхностная плотность заряда. Полученные решения приведены в сравнении с результатами других авторов.

3.2.4. Контактная задача Герца для электроупругого полупространства [130] Рассмотрим гранично-элементные решения статической начально-краевой задачи о действии контактной нагрузки Герца F(r) Q(r ),0, P(r ), где P( r ) p0 1 r 2 c 2, p0 121 109 Па, Q(r ) fP (r ), f 0.25, r x12 x2, на часть дневной поверхности электроупругого полупространства (рис. 176). Дневная поверхность является свободной – на ней заданы поверхностные усилия ti 0 Па (i 1,3) и поверхностная плотность заряда q 0 Кл м. Нагружаемый участок ограничен окружностью радиусом c 1 м с центром в начале координат.

В качестве материала взят трансверсально изотропный цирконат-титанат свинца – пьезокерамик PZT5A со следующими параметрами [130]:

Гранично-элементная сетка строится с учетом одной плоскости симметрии.

Половина расчетной сетки содержит 696 элементов и 741 узел (рис. 177).

На рис. 178, 179 представлены полученные по формуле Сомильяны величины откликов нормализованных напряжений t3 p0 и поверхностной плотности заряда qc11 e33 p0 для точек с координатами x1 c [4,1] [1,4], x2 c 0, x3 c 2.

t3/p qc11/e33p полупространство ослаблено полостью, расположенной под площадкой нагружения (рис. 280). Форма полости определяется следующим уравнением:

Гранично-элементная сетка строится с учетом одной плоскости симметрии.

Половина расчетной сетки содержит 1128 элементов (рис. 181).