Пособие учебно-методическое Рекомендовано методической комиссией механико-математического факультета для студентов ннгу, обучающихся по направлению подготовки 010500 «Прикладная математика и информатика» МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Нижегородский государственный университет им. Н.И. Лобачевского С.А. Капустин Метод взвешенных невязок решения задач механики деформируемых тел и теплопроводности Учебно-методическое пособие Рекомендовано методической комиссией механико-математического факультета для студентов ННГУ, обучающихся по направлению подготовки 010500 «Прикладная математика и информатика» Нижний Новгород 2010 УДК 539.3 (075) ББК В25 (Я73-4) К20 К 20 Капустин С.А. МЕТОД ВЗВЕШЕННЫХ НЕВЯЗОК РЕШЕНИЯ ЗАДАЧ МЕХАНИКИ ДЕФОРМИРУЕМЫХ ТЕЛ И ТЕПЛОПРОВОДНОСТИ: учебно-методическое пособие. – Нижний Новгород: Нижегородский госуниверситет, 2010. – 60 с. Учебно-методическое пособие посвящено изложению основ наиболее известных методов дискретизации непрерывных задач, определенных соответствующими дифференциальными уравнениями, которые являются основой для построения широкого класса методов численного решения задач математической физики. В частности, рассмотрены метод конечных разностей и метод взвешенных невязок, объединяющий в своем составе ряд различных конкретных методов, таких как метод коллокации, метод Галеркина, метод конечных элементов и др. Пособие рассчитано на студентов бакалавриата и магистратуры университетов по специальностям «Прикладная математика» и «Механика». УДК 539.3 (075) ББК В25 (Я73-4) © Нижегородский государственный Университет им. Н.И. Лобачевского, 2010 Содержание Введение. 4 1. Общая характеристика уравнений теории упругости и теплопроводности. Метод конечных разностей. 7 1.1. Уравнения теории упругости. 7 1.2. Уравнения теплопроводности. 10 1.3. Конечно-разностные аппроксимации производных. 12 1.4. Решение одномерных задач методом конечных разностей. 14 2. Метод взвешенных невязок с использованием базисных функций. 20 2.1. Аппроксимация функций с использованием систем базисных функций. 20 2.2. Основы метода взвешенных невязок. 22 2.3. Аппроксимация решений дифференциальных уравнений. 25 2.4. Использование в МВН функций, не удовлетворяющих априори краевым условиям.. 28 2.5. Естественные краевые условия. 31 2.6. Общая формулировка естественных краевых условий для задач теплопроводности. 33 2.7. Применение МВН в задачах теории упругости. 34 3. Метод взвешенных невязок с кусочным определением базисных функций (метод конечных элементов) 37 3.1. Особенности задания базисных функций в методе конечных элементов (МКЭ) 37 3.2. Аппроксимация решения дифференциальных уравнений с использованием кусочно-определенных базисных функций. 38 4. Построение базисных координатных функций в МКЭ.. 45 4.1. Основные требования к координатным функциям в МКЭ.. 45 4.2. Построение базисных функций конечных элементов в обобщенных координатах. 45 4.3. Построение локальных систем координат. 48 4.4. Лагранжево семейство элементов. 50 4.5. Сирендипово семейство элементов. 52 4.6. Треугольные элементы.. 53 Литература. 60 Введение Уровень развития современной техники неразрывно связан с состоянием исследований в области механики и вычислительных методов, позволяющих осуществлять математическое моделирование реальных физических процессов. Поэтому количественные методы исследования таких процессов проникают в настоящее время практически во все сферы человеческой деятельности, а математические модели становятся основным средством познания реальной действительности. Количественное описание исследуемых процессов с использованием математических моделей обычно строится путем последовательной реализации двух основных этапов. На первом этапе осуществляется математическая постановка решаемой задачи, заключающаяся в формулировке исходных дифференциальных уравнений, ограниченных определенной системой краевых и начальных условий. Далее требуется решить поставленную задачу для конкретных конфигураций исследуемых объектов, заданных свойств этих объектов, конкретных начальных и краевых условий. Решение таких задач на основе современных аналитических методов обычно удается получить лишь для простейших систем уравнений, ограниченных видов конфигураций объектов и частных видов краевых условий. Кроме этого, применение аналитических методов не позволяет в должной мере использовать для решения задач постоянно растущие возможности современных ЭВМ. Более перспективными в этом плане являются методы вычислительной математики, непосредственно ориентированные на применение ЭВМ и постоянно совершенствующиеся вместе с прогрессом в области вычислительной техники. На основе этих методов удается преобразовать исходную задачу к чисто алгебраической форме, для решения которой достаточно использовать только основные арифметические операции. Такие преобразования возможны на основе применения различных методов дискретизации непрерывных задач, определенных соответствующими дифференциальными уравнениями. В результате дискретизации бесконечные множества чисел, определяющих неизвестные функции, аппроксимируются конечным числом параметров, а исходная непрерывная задача сводится к определению значений этих параметров. К настоящему времени известно большое число различных методов дискретизации, отличающихся различными способами аппроксимации искомых функций и исходной непрерывной задачи. Применительно к решению краевых задач математической физики, вопросам решения которых посвящен основной материал предлагаемого пособия, в настоящее время наибольшее распространение получили две наиболее известные группы методов дискретизации, существенно различных по своей основе, но, тем не менее, имеющих много общих черт. Первую группу составляют конечно-разностные методы (МКР), имеющие в зависимости от постановки решаемых задач большое число различных модификаций. В самых общих чертах суть МКР сводится к замене дифференциальных выражений в исходных дифференциальных уравнениях некоторой краевой задачи соответствующими конечно-разностными выражениями относительно значений искомых функций в узлах сеточной области. С этой целью исследуемая область покрывается системой сеточных линий. Действительная граница области заменяется сеточной границей таким образом, чтобы сеточные узлы наилучшим образом приближали истинную границу. В каждом узле области дифференциальные операторы краевой задачи заменяются соответствующими конечно-разностными операторами. При этом исходная краевая задача сводится к задаче решения системы алгебраических уравнений относительно значений искомых функций в узлах сеточной области. Вторая группа методов с общим названием “метод взвешенных невязок” (МВН) включает в себя широкий класс различных схем, в основу которых положено два, общих для всех конкретных схем, фактора: - введение системы базисных функций, удовлетворяющих определенным условиям исходной задачи и используемых для аппроксимации искомых функций; - построение метода определения параметров аппроксимации путем приравнивания нулю интеграла от невязки исходных уравнений, взятых с определенной системой весовых функций. Выбирая различные варианты базисных функций, весовых функций, различные способы построения уравнений невязок для определения параметров аппроксимации искомых функций можно, получать различные варианты конкретных схем МВН. Обычно базисные функции выбирают такими, чтобы для них выполнялись те или иные условия исходной краевой задачи, причем, для определения параметров этих функций, составляются уравнения невязок для оставшейся части невыполненных условий. В частности, при использовании базисных функций, удовлетворяющих всем краевым условиям исходной задачи, для определения параметров этих функций составляется уравнение в виде равенства нулю интеграла по области от произведения невязки основного дифференциального оператора задачи на выбранную систему весовых функций. В результате в зависимости от выбора весовых функций могут быть получены соответствующие классические схемы МВН: метод моментов, методы коллокации, метод Галеркина и другие. При выборе базисных функций, удовлетворяющих основному дифференциальному оператору внутри области, для определения неизвестных параметров составляется уравнение, содержащее интеграл по границе области от произведения невязки граничных условий на соответствующую систему весовых функций. Получаемые при этом конкретные схемы МВН соответствуют схеме метода граничных уравнений. В классических вариантах МВН базисные функции задаются для всей области, что приводит к определенным трудностям при выборе самих функций и при численном решении таких задач. В частности известно, что с увеличением номера приближения в таких схемах существенно ухудшается обусловленность результирующей системы линейных алгебраических уравнений. Однако перечисленные трудности могут быть легко преодолимы путем использования локально определенных базисных функций, задаваемых на отдельных подобластях исследуемой области. Применение такой системы МВН, которую в дальнейшем будем называть методом конечных элементов (МКЭ), предполагает деление исследуемой области системой непересекающихся подобластей (конечных элементов), в каждой из которых осуществляется аппроксимация искомых функций с помощью простейших степенных полиномов. В результате существенно упрощается сам выбор базисных функций, удовлетворяющих заданным краевым условиям, появляются широкие возможности исследования областей, имеющих нерегулярные границы, улучшается структура и обусловленность результирующей системы. Следует иметь в виду, что для задач, допускающих вариационную формулировку исходных уравнений, идентичные МВН схемы МКЭ могут быть получены в рамках известных вариационных методов. В частности, в первых работах, определивших появление МКЭ, этот метод трактовался либо как разновидность классического вариационного метода Релея-Ритца с локальным заданием базисных функций, либо как разновидность вариационно-разностного метода. Поэтому МКЭ можно рассматривать как общий метод дискретизации непрерывных краевых задач, имеющий непосредственную связь с МВН и классическими вариационными методами, использующими базисные функции для аппроксимации неизвестных функций, а также с конечно-разностными методами, основанными на вариационной формулировке исходных уравнений. Вобрав в себя лучшие качества конечно-разностных методов, МКЭ имеет перед ними определенные преимущества, обеспечивая возможность свободного размещения узловых точек и сгущения их в местах высокого градиента функций, возможность одновременного использования в рамках единой расчетной схемы конечных элементов различной сложности и мерности. Кроме того, наличие разрывов в геометрии конструкции нисколько не затрудняет анализ таких конструкций МКЭ, что делает этот метод весьма перспективным для исследования сложных конструктивных образований. В настоящее время МКЭ признается одним из наиболее универсальных методов решения прикладных задач математической физики, успешно используемых для математического моделирования в области прочности, аэрогидродинамики, теплофизики. Предлагаемое учебное пособие посвящено изложению основ перечисленных методов дискретизации и иллюстрации их применения для двух классов задач: -задач стационарной теплопроводности, на примерах решения которых удается наиболее просто и наглядно продемонстрировать основные принципы построения численных решений; -задач теории упругости, являющихся наиболее сложными задачами математической физики, и требующих для построения численных решений привлечение более широкого набора вычислительных средств. В первой главе пособия рассмотрены постановка задач и основные уравнения теории упругости. Приведена система дифференциальных уравнений в частных производных и краевых условий для функций перемещений, описывающих упругое равновесие деформируемых систем. Рассмотрены основные соотношения и совокупность уравнений, описывающих распределение температур в изотропных теплопроводящих средах при различных граничных условиях. Дано описание способов построения конечно-разностных аппроксимаций производных различных порядков и рассмотрены вопросы применения МКР для решения одномерных задач. Вторая глава посвящена изложению основ МВН с использованием базисных функций, определенных на всей области рассматриваемой задачи. На примерах решения задач теплопроводности рассмотрены различные варианты аппроксимации дифференциальных уравнений на основе различных систем весовых функций и базисных функций, удовлетворяющих и неудовлетворяющих априори заданным краевым условиям. Рассмотрена слабая формулировка МВН для задач теплопроводности и теории упругости. Третья глава содержит описание основ МКЭ, как специальной формы МВН с кусочным определением базисных функций. Рассмотрены особенности задания локальных базисных функций в МКЭ и применения их для аппроксимации дифференциальных уравнений задач теплопроводности и теории упругости. Четвертая глава посвящена вопросам построения базисных координатных функций в МКЭ. Сформулированы основные требования к таким функциям, обеспечивающие сходимость численных решений. Рассмотрены некоторые способы построения базисных функций в элементах различной конфигурации и пространственной мерности. Приведены основные сведения об искривленных изопараметрических элементах и применении для вычисления их жесткостных характеристик схем численного интегрирования. 1. Общая характеристика уравнений теории упругости и теплопроводности. Метод конечных разностей 1.1. Уравнения теории упругости Математическая теория упругости изучает количественные соотношения, характеризующие деформации или внутренние относительные смещения в твердом теле при заданных внешних воздействиях, в виде объемных сил, распределенных по объему тела V, а также поверхностных сил и граничных перемещений, заданных на поверхностях и , соответственно. В теории упругости тело рассматривается как непрерывная сплошная среда, для которой считаются справедливыми гипотезы упругости (способности тела полностью восстанавливать свою форму после устранения причин, вызвавших деформацию), однородности (независимости свойств материала в пределах рассматриваемого объема тела) и линейной зависимости между возникающими в теле напряжениями и деформациями. Пусть некоторое тело, занимающее объем V , ограниченный поверхностью , находится в условиях статического равновесия при постоянной температуре Т= const . под действием объемных сил , поверхностных сил , а также граничных перемещений , заданных на части поверхности . В условиях перечисленных воздействий точки тела получают перемещения , удовлетворяющие условию (1.1) Кроме того, за счет смещения отдельных точек относительно друг друга, в теле возникнут деформации, характеризуемые симметричным тензором 2ого ранга и связанные с перемещениями точек тела известными дифференциальными соотношениями. В случае если перемещения и деформации в теле малы, то эти соотношения линейны и имеют вид . (1.2) В литературе зависимости (1.2) известны как соотношения Коши. В матричной форме эти соотношения можно представить в виде , (1.3) где - векторы компонент тензора деформаций (в силу симметрии тензора в вектор включены лишь 6 независимых компонентов, причем сдвиговые компоненты приняты в виде ) и перемещений в точке (1.4) [d ] - матричный дифференциальный оператор . (1.5) Появление деформаций в теле в свою очередь вызывает появление в нем силовых полей, характеризуемых тензорным полем напряжений . В упругом теле напряжения связаны с деформациями известными соотношениями упругости (обобщенный закон Гука): . (1.6) В общем случае число независимых компонент тензора (с учетом симметрии тензоров ) равно 21. В случае изотропных тел число независимых величин, определяющих компоненты , сокращается до 2 , (1.7) где e - шаровая составляющая тензора деформации , - упругие постоянные Ламе ( - модуль деформации сдвига) , (1.8) где - модуль объемной деформации, - коэффициент поперечной деформации. В матричной форме соотношения (1.6) могут быть записаны в виде , (1.9) где - вектор компонент тензора напряжений (содержащий 6 независимых компонент) . (1.10) - матрица коэффициентов упругости материала, составленная из коэффициентов . Для упругого изотропного материала эта матрица имеет вид (1.11) Напряжения в теле, находящемся в равновесии под действием сил , и граничных смещений , должны удовлетворять дифференциальным уравнениям равновесия в объеме тела V : (1.12) и уравнениям равновесия на части границы , , (1.13) где - направляющие косинусы нормали к поверхности границы . В матричной форме эти уравнения могут быть переписаны в виде , (1.14) , (1.15) где - матричный дифференциальный оператор, определенный соотношением (1.5), {F } и {P } - векторы объемных и поверхностных сил (1.16) - матрица направляющих косинусов нормали к поверхности (1.17) Приведенные уравнения (кинематические (1.1) − (1.3), физические (1.6) − (1.11) и статические (1.12) − (1.17)) составляют полный комплект уравнений, статических задач линейной теории упругости. При непосредственном решении таких задач приведенные уравнения обычно преобразуются к некоторому виду разрешающих уравнений, в которых в качестве основных неизвестных фигурируют поля какого-либо одного типа, например, напряжений или перемещений . В частности для решения задач в перемещениях в уравнениях (1.12) и (1.13) напряжения выражаются через перемещения с помощью соотношений (1.2) и (1.6). (1.18) где . (1.19) В результате решение задачи теории упругости сводится к решению системы дифференциальных уравнений в частных производных 2-го порядка относительно функций перемещений при граничных условиях, заданных на частях поверхности и (1.20) , (1.21) где - единичная матрица. 1.2. Уравнения теплопроводности Рассмотрим закономерности распределения температуры T в теле, занимающем объем V , ограниченный поверхностью . Баланс тепла в единичном объеме за единицу времени соответствует дифференциальному уравнению , (i=1-3 ), (1.22) где qi – поток тепла в направлении x i на единицу длины за единицу времени, Q= Q (xi , t ) – тепло, генерируемое в элементарном объеме за единицу времени (источник тепла), - изменение тепла за счет изменения температуры T =T (xi , t ), , c – плотность и теплоемкость материала. В свою очередь поток тепла в изотропной среде в произвольном направлении n связан с изменением температуры T соотношением , (1.23) где k – коэффициент теплопроводности материала. Записывая выражения для потоков тепла в направлениях xi ( ) можно привести уравнение баланса (1.22) к виду . (1.24) Уравнение (1.24) описывает задачи распределения температур в пространственно-временной области, характеризуемой независимыми переменными xi , t . Для решения задач нестационарной теплопроводности к уравнению (1.24) необходимо добавить начальные условия, характеризующие распределение температуры в некоторый начальный момент времени и граничные условия на границе области , которые сводятся к двум основным типам: - заданию температуры на части границы (1.25) - заданию потока на части в направлении нормали к границе n . (1.26) В случае установившихся процессов теплопроводности распределение температур перестает зависеть от времени и уравнение (1.24) упрощается , (1.27) где T= T( xi ), Q= Q( xi ). Решение уравнения (1.27) стационарной теплопроводности требует задания лишь граничных условий типа (1.25) и (1.26). Формальную запись совокупности уравнений, описывающих стационарную задачу теплопроводности, можно представить в виде: В общем случае задача теплопроводности (1.25) - (1.27) является нелинейной, так как коэффициент теплопроводности k может зависеть от температуры. Если такая зависимость отсутствует и k=const , то уравнение (1.27) становится линейным . (1.28) При исследовании процессов теплопроводности в двумерной постановке уравнения (1.25), (1.26), (1.28) приобретают вид (V – двумерная область, занимаемая телом), (1.29) ( - части контура области V ). В одномерном случае аналогичные уравнения сводятся к выражениям (V –отрезок оси x , L – длина отрезка) и (или) (1.30) 1.3. Конечно-разностные аппроксимации производных Среди численных методов решения задач математической физики широкое применение находят разностные методы. При использовании этих методов рассматриваемая область аппроксимируется сеточной областью, а значения производных искомых функций заменяются разностными отношениями через значения этих функций в узлах сетки. Далее для каждого узла сеточной области составляются соответствующие разностные аналоги исходных функциональных уравнений, составляются разностные аналоги заданных граничных условий и задача сводится к решению системы алгебраических уравнений. Для построения простейших разностных аппроксимаций функций в одномерном случае запишем выражение для функции f (x ) в малой окрестности x +Δx с помощью разложения ее в ряд Тейлора (1.31) Пусть функция f (x ) задана в виде таблицы , . Тогда, ограничиваясь двумя членами разложения, можно записать , где или (1.32) Погрешность этой формулы имеет порядок . Такое представление производной носит название аппроксимации первой производной разностью вперед. Аналогичное выражение можно получить для аппроксимации производной разностью назад (1.33) Для получения более точных формул можно представить выражение для функций в узлах i -1 и i +1 в виде , (1.34) , (1.35) Вычитая эти равенства одно из другого можно получить . (1.36) Погрешность E этой аппроксимации определяется . Удерживая в выражениях (1.34), (1.35) слагаемые, содержащие и складывая результаты можно получить разностное выражение для второй производной , (1.37) где . Подобным образом можно построить разностные выражения для производных высших порядков ; (1.38) . (1.39) Аналогичная техника может быть использована для построения разностных аппроксимаций частных производных в двумерном и пространственном случаях. В частности на двумерной сетке в окрестности узла с индексами i, k частные производные второго порядка функции f (x, y ) могут быть записаны в виде: ; (1.40) ; (1.41) . (1.42) 1.4 . Решение одномерных задач методом конечных разностей В качестве примера решения методом конечных разностей (МКР) одномерных уравнений второго порядка рассмотрим решение задачи теплопроводности в одномерном случае. Пусть требуется определить функцию T (x ), удовлетворяющую дифференциальному уравнению на отрезке 0< x< L . Граничные условия задачи соответствуют заданию на каждом конце отрезка температуры или потока . Для решения задачи МКР отрезок, на котором ищется функция, разбивается определенным числом равностоящих точек с координатами . При этом . Далее для всех внутренних точек сетки составляется конечно-разностный аналог исходного дифференциального уравнения на основе применения полученных выше формул конечно-разностного представления второй производной. В типичном узле это уравнение будет иметь вид . (1.43) Если на обоих краях отрезка заданы значения температуры , то в первом и последнем узлах уравнения не составляются, а во втором и предпоследнем уравнениях известные значения температуры переносятся в правую часть. В результате получается система L -1 алгебраических уравнений (АУ) относительно L -1 неизвестных значений в узлах сеточной области (1.44) В матричном представлении полученную систему уравнений можно записать в виде , (1.45) где - вектор неизвестных узловых значений температур , [K ] – матрица системы (1.46) {R } – вектор правых частей системы (1.47) Следует отметить, что матрица [K ] получается симметричной и узколенточной (трехдиагональной). Таким образом, исходная задача определения неизвестных непрерывных функций заменяется задачей решения системы АУ относительно дискретных значений T 1 . . .TL -1 . Иначе говоря, МКР дает информацию о значениях функций в узлах сеточной области, но не дает информации о значениях функций между точками, т.е. дифференциальное уравнение аппроксимируется только в конечном числе точек. Рассмотрим также случай, при котором на одном из концов отрезка (например, при x =L ) задан поток тепла при x =L. (1.48) При этом, поскольку значение температуры TL оказывается неизвестным для определения всех узловых температур к системе (1.44) необходимо добавить еще одно уравнение. Такое уравнение можно получить, записав в конечных разностях уравнение (1.48) в узле x =L . Если аппроксимировать производную разностью назад, то уравнение (1.48) примет вид (1.49) Добавляя уравнение (1.49) к системе (1.44) получим систему L уравнений, необходимых для определения L неизвестных значений температур. (1.50) Полученная таким образом система обладает существенным недостатком, связанным с различным порядком аппроксимации решения внутри области и на границе . Этого недостатка можно избежать, если добавить к системе (1.44) уравнение (1.43) в узле x =L и использовать для условия (1.48) аппроксимацию производной в центральных разностях. В результате система алгебраических уравнений примет вид (1.51) где TL +1 – значение температуры в законтурной точке. В качестве примера решения МКР уравнений четвертого порядка рассмотрим задачу изгиба балки, нагруженной равномерно распределенной по длине поперечной нагрузкой. Дифференциальное уравнение изгиба балки , (1.52) где - функция поперечного перемещения (прогиба) балки, p - интенсивность поперечной нагрузки, D – жесткость балки на изгиб, l - длина балки. В качестве граничных условий задачи на каждом из торцов балки могут быть заданы следующие варианты условий: - перемещение W балки на опоре или значение перерезывающей силы ; - угол поворота или значение изгибающего момента на торце балки . По аналогии с задачей теплопроводности, для решения задачи отрезок оси x, на котором ищется функция W (x ) разбивается системой равноотстоящих точек на n отрезков с координатами . Кроме этого для реализации граничных условий на концах балки добавляются еще по одной или по две законтурные точки . Для каждого из внутренних узлов сетки составляется конечно-разностный аналог дифференциального уравнения (1.52) с использованием конечно-разностного представления производной четвертого порядка. В типовом узле i такое уравнение будет иметь вид ; (1.53) . Аналогичные уравнения записываются и для граничных узлов, если в них значения функций W являются неизвестными. Уравнения, содержащие известные значения прогибов в граничных узлах, корректируются путем перенесения известных слагаемых в правую часть. Далее осуществляется конечно-разностная аппроксимация граничных условий, записанных в виде производных для функции W. Получающиеся при этом разностные уравнения позволяют определить значения функции W в законтурных узлах, появляющихся при записи уравнений в узлах . Примеры Задача 1 Найти решение уравнения (0< x< 1.) при краевых условиях x =0, T0 =0; x =1, TL =1. Используем для решения сетку n =3; шаг сетки Согласно краевым условиям T0 =0., T3 =1. Неизвестны T1 и T2 Уравнение в точке x= xi имеет вид Конечно-разностный аналог уравнения в узле i: Ti +1 -2Ti +Ti -1 -Ti Δx 2 =0 Система АУ принимает вид (уравнения составляются в узлах i =1, 2) -1/9T 1 +T 2 =0 T 1 -2T 2 +T 3 -1/9T 2 =0 или Решение системы: T 1 =0.2893; T 2 =0.6107 Точное аналитическое решение T 1 =0.2889; T 2 =0.6102. Задача 2 Найти решение уравнения задачи 1 при краевых условиях x =0; T =0. x =1; . При использовании рассмотренной выше конечно-разностной сетки в задаче известны T0 =0, неизвестны T1 , T2 , T3 . Конечно-разностные уравнения для узлов i =1 и 2 остаются без изменения. Для узла i =3 составим уравнение для потока с использованием разности назад (погрешность аппроксимации. 0(Δx )): (T2 -T3 )/(1/3)=1. Результирующая система АУ примет вид Точное аналитическое решение T 1 =0.220, T2 =0.4648, T 3 =0.7616. Задача 3 Найти решение задачи 1 при краевых условиях задачи 2 с использованием конечно-разностной аппроксимации с погрешностью 0(Δ x )2 . Первые два конечно-разностных уравнения остаются такими же, как в примере 2. Добавляются уравнения в узле i =3: =0 (T 4 – значение температуры в законтурной точке) и уравнение для потока в узле i =3 с использованием производной в центральных разностях (погрешность аппроксимации 0(Δx )2 ) Результирующая система АУ принимает вид Результаты получаются более точными, но нарушается симметрия матрицы. Задача 4 Вычислить функцию прогиба W( x) балки длиной l =1 м., жестко заделанной по торцу x=0 и опертой на непроседающую опору по торцу x= l , находящейся под действием равномерно распределенной поперечной нагрузки p и опорного момента . Жесткость балки на изгиб ; p =1 н/м ; Дифференциальное уравнение изгиба балки Граничные условия Схема дискретизации Конечно-разностный аналог уравнения изгиба для внутреннего (i- го) узла можно представить в виде ; Конечно-разностные уравнения на левой опоре x =0: в узле “0” уравнение не составляется . Конечно-разностные уравнения на правой опоре x= l : ; в узле n уравнение не составляется ; , . или . Результирующая система имеет вид Неизвестные функции Общее число неизвестных 2. Метод взвешенных невязок с использованием базисных функций 2.1. Аппроксимация функций с использованием систем базисных функций Рассмотренный выше метод конечных разностей позволяет свести решение дифференциальных уравнений для непрерывных функций к численному решению системы алгебраических уравнений путем разностной аппроксимации производных в узлах сеточной области. В результате решения такой задачи удается определить значения искомой функции в конечном числе точек. Однако существует много других способов аппроксимации функций, которые могут быть использованы для численного решения дифференциальных уравнений. В частности широкое распространение получили методы, основанные на применении базисных функций, определенных на всей области решения поставленной задачи [1]. Пусть требуется аппроксимировать заданную функцию f в области V , ограниченной поверхностью Г и принимающую на этой поверхности заданные значения . С этой целью введем функцию , удовлетворяющую условию , а также систему линейно-независимых базисных функций , удовлетворяющих на граничной поверхности Г условиям для всех m . Тогда аппроксимацию функции f можно представить в виде , (2.1) где (m =1,2…M ) - некоторые константы, определяемые из условия наилучшего приближения функции . Система базисных функций должна обеспечивать улучшение аппроксимации функции (т.е. сходимость ) при увеличении числа М базисных функций. Для этого комбинация (1.1) должна удовлетворять условиям полноты, т.е. обеспечивать возможность сколь угодно точного представления произвольной функции f , удовлетворяющей условиям на границе . Таким образом, для получения хорошей аппроксимации функции с помощью комбинации типа (1.1) необходимо решить следующие вопросы: - подобрать функцию , удовлетворяющую условию ; - выбрать систему базисных функций Nm , удовлетворяющих условиям полноты; - выбрать способ определения констант при выбранных базисных функциях Nm . Первый вопрос обычно решается сравнительно просто, поэтому ниже основное внимание будет направлено на вопросы выбора базисных функций и способа вычисления констант. В настоящее время известен ряд классов функций, которые могут быть использованы в качестве базисных в аппроксимациях типа (2.1). В частности широко используются базисные функции в виде степенных полиномов, приводящие к аппроксимации функций алгебраическими многочленами вида , (2.2) ; . Многочлены типа (2.2) хорошо изучены и удобны в использовании: легко вычисляются, без труда дифференцируются и интегрируются. Для вычисления параметров может быть использован способ интерполяции, согласно которому искомая функция заменяется аппроксимирующей, совпадающей с искомой в ряде различных точек (узлах интерполяции). При этом параметры интерполяции определяются из решения системы алгебраических уравнений ; ; ; или . (2.3) Эта система однозначно разрешима, если узлы интерполяции попарно различны. Далее аппроксимацию (2.2) можно привести к виду ; ; . Следует отметить, однако, что на практике такой способ определения параметров и представления аппроксимации вида (2.2) используется достаточно редко, т.к. функции Nm определены неявно из решения системы (2.3), которая обычно оказывается плохо обусловленной. В связи с этим более широкое применение находят другие формы аппроксимации функций алгебраическими многочленами, в которых функции Nm определены в явном виде, например в форме многочленов Лагранжа , (2.4) где . Особенностью аппроксимации типа (2.4) является то, что в ней принимается , а базисные функции обладают свойством (2.5) В качестве базисных функций также часто используются тригонометрические функции, приводящие к аппроксимации вида , (2.6) где - длина отрезка , - линейная функция, принимающая на границе отрезка значения . На основе теории рядов Фурье коэффициенты этой аппроксимации могут быть получены в виде . (2.7) Возможны и другие частные способы построения аппроксимаций типа (1.1). При этом существует общий метод определения констант в аппроксимации (1.1), на основе которого можно получить все рассмотренные выше подходы. Этот метод носит название – метод взвешенных невязок [1]. 2.2. Основы метода взвешенных невязок Наиболее общий метод определения параметров в аппроксимации вида (2.1) можно получить, если ввести невязку аппроксимации Rv (2.8) Поскольку Rv представляет собой функцию, зависящую от координат точек области V , для уменьшения величины этой невязки, т.е. для приближения функции потребуем выполнения условий , (2.9) где {Wn ; n =1,2,…M } – множество линейно-независимых весовых функций. Тогда, при условии полноты выбранных весовых функций Wn , сходимость функций будет достигнута, если потребовать выполнение равенств (2.9) для всех n при . Подставляя в (2.9) представление функции согласно (2.1) можно получить или (2.10) Уравнение (2.10) представляет собой систему линейных алгебраических уравнений, из решения которой могут быть найдены неизвестные константы . По аналогии с предыдущим разделом эту систему можно записать в матричном виде где , (2.11) , . Задаваясь различными видами весовых функций на основе (2.10) можно получать различные варианты метода взвешенных невязок. Ниже рассмотрены некоторые наиболее часто используемые формы представления весовых функций [1]. Коллокация в точке. В этом случае полагается, что , где - дельта функция Дирака, обладающая свойствами , . При этом . Иначе говоря, принимается, что Wn =1 в точке xn и Wn =0 нулю во всех других точках. Согласно такому выбору весовой функции невязка Rv оказывается равной нулю в ряде заданных точек xn , а элементы системы (2.11) принимают вид . Понятно, что на основе схемы поточечной коллокации могут быть получены различные способы интерполяции функций. Коллокации по подобластям. В этом методе принимается, что Wn =1 в некоторой подобласти и W =0 при . Иначе говоря, принимается условие, согласно которому интеграл от невязки обращается в ноль для ряда подобластей . Элементы системы (2.11) при этом принимают вид ; Метод Галеркина. В качестве весовых функций выбираются базисные функции Wn =Nn . Элементы системы уравнений (2.11) при этом приобретают вид . Особенностью метода Галеркина является симметричность матрицы [K ]. Кроме этого, если в качестве Nm использовать систему ортогональных функций, таких, что , то матрица [K ] становится диагональной. Например, если на отрезке 0<x <L использовать для аппроксимации систему базисных функций , то ; . При этом могут быть сразу найдены коэффициенты : . Такая форма аппроксимации соответствует приближению функций на основе рассмотренных выше рядов Фурье. Метод моментов. Весовые функции принимаются в виде Рис. 2.1. Весовые функции для момента нулевого порядка - момент нулевого порядка – суммарная площадь(см.рис.2.1) - момент первого порядка – момент площадей относительно начала координат; - момент второго порядка и т.д. Метод наименьших квадратов. В методе наименьших квадратов минимизируется интеграл от квадрата невязки функций по области V . . Для нахождения минимума I используется условие стационарности, приводящее к системе уравнений . Принимая для функции представление в виде (2.1) и учитывая, что , условие стационарности функционала I можно привести к виду Полученные уравнения совпадают со стандартной формой метода взвешенных невязок с весами Wn = Nn .В рассматриваемом случае формулировка метода наименьших квадратов совпадает с методом Галеркина. 2.3. Аппроксимация решений дифференциальных уравнений Рассмотренный выше общий метод аппроксимации функций может быть легко распространен на аппроксимацию решений дифференциальных уравнений. Рассмотрим дифференциальное уравнение в области V для некоторой функции , (2.12) где D – линейный дифференциальный оператор, P =P (xi ) – заданная функция координат области V . Решение уравнения (2.12) должно удовлетворять краевым условиям на границе области Г , которые в общем виде могут быть записаны , (2.13) где G – линейный оператор, а р (Г ) – заданная функция координат границы Г . Например, для рассмотренных выше двумерных задач теплопроводности функции , , а уравнения (2.12) и (2.13) примут вид ; (2.14) (2.15) По аналогии с предыдущим разделом искомую функцию f можно аппроксимировать с помощью базисных функций , (2.16) удовлетворяющих краевым условиям на границе области Г . (2.17) Если функции Nm непрерывны в V и все их производные существуют, то на основе (2.16) можно получить аппроксимации производных функций f ; (2.18) . Поскольку представление функции f в виде (2.16) удовлетворяет краевым условиям на Г , для получения аппроксимации f вычислим невязку Rv уравнения (2.12) . (2.19) Выбирая по аналогии с предыдущим разделом систему весовых функций {Wn ; n =1,2,..} потребуем выполнения условия Rv =0 в V в виде . (2.20) Поскольку Wn , Ψ , Nm и P являются заданными функциями координат, систему уравнений (2.20) можно привести к системе алгебраических уравнений ; ; (2.21) ; . Если каждая из функций Wn и Nm определена на всем пространстве области V , то в общем случае матрица системы (2.21) оказывается заполненной. Для получения конкретных значений элементов системы (2.21) необходимо выбрать соответствующие системы базисных Nm и весовых Wn функций, причем от того, насколько удачно выбраны эти функции, будет зависеть качество и общая эффективность численного решения. Пример. Рассмотрим пример решения задачи теплопроводности f =T (x ) для отрезка , при краевых условиях T=0 при x=0 и T=1 при x=1, или Gf= T, =0 при x =0 и p = -1 при x =1. В качестве функции можно выбрать функцию , удовлетворяющую краевым условиям при x =0 и x =1. В качестве базисных функций выберем систему , удовлетворяющую условиям Nm =0 при x =0 и x =1 для всех m . Тогда предложенная выше комбинация будет удовлетворять заданным краевым условиям и может быть использована для аппроксимации искомой функции. Рассмотрим два варианта выбора весовых функций Wn при M =2: - поточечную коллокацию при x 1 =1/3 и x 2 =2/3; - метод Галеркина ; Метод поточечной коллокации. ; ; n =1 m =1 ; m =2 ; R1 =1/3; n =2 m =1 ; m =2 . В результате решения полученной системы могут быть найдены параметры а m a1 = -0.05312; a2 = 0.004754. Подставляя полученные значения параметров в приближенное представление функции Т можно вычислить значения температуры в узлах x =1/3 и x =2/3: T1 = 0,2914 и T2 = 0,6165. Метод Галеркина. n =1 m =1 ; m =2 ; R1 = ; n =2 m =1 ; m =2 R2 = . Из решения полученной системы можно получить a1 = -0,05857; a2 = 0,007864 Значения температуры при этом в узлах x 1 =1/3 и x 2 =2/3 получаются равными Т 1 = 0,2894; T 2 = 0,6091. Для сравнения в таблице 2.1 приведены результаты решения этой задачи, полученные на основе поточечной коллокации, метода Галеркина (МГ), метода конечных разностей (МКР) и точного решения (ТР) для двух точек отрезка x =1/3 и x =2/3. Таблица 2.1 Значения температур для х =1/3 и х =2/3 на основе точного решения (ТР), методов коллокации (ПК), Галеркина (МГ) и конечных разностей (МКР) X ПК МГ МКР ТР 1/3 0,2914 0,2894 0,2893 0,2889 2/3 0,6165 0,6091 0,6107 0,6102 Метод наименьших квадратов ; Таким образом, в методе наименьших квадратов Wn = DNn - отличается от метода Галеркина, в котором Wn = N n . 2.4. Использование в МВН функций, не удовлетворяющих априори краевым условиям В связи с трудностью подбора аппроксимирующих функций, удовлетворяющих всем краевым условиям, можно потребовать выполнения этих условий приближенно, в рамках метода взвешенных невязок. В частности аппроксимирующую функцию f можно представить в виде , (2.22) где не удовлетворяет заранее каким-либо краевым условиям. Для приближенного удовлетворения этих условий к невязке Rv дифференциального уравнения можно добавить невязку краевых условий RГ (2.23) и потребовать выполнения условий (2.24) где - некоторая система весовых функций, выбираемая по аналогии с Wn . Подставляя (2.22) и (2.23) в (2.24) можно получить систему алгебраических уравнений типа (2.21) для определения параметров . ; ; (2.25) . Например, для решения рассмотренной выше задачи теплопроводности с краевыми условиями Т =0 при x =0 и T =1 при x =1 можно использовать систему базисных функций . В данном случае граничные условия задаются для двух точек x =0 и x =1. При этом условие (2.24) приводится к виду . Если в качестве весовых функций принять , то уравнение МВН для рассматриваемой задачи может быть записано ; ; . При M =1; ; ; R 1 =-1; =1/3; При M =3; ; N 1 =1; N 2 =x ; N 3 =x 2 ; ; ; ; и т.д. ; ; . В результате решения системы a1 = 0,068; a2 = 0,632; a3 = 0,226 или . Для сравнения в таблице 2.2 приведены результаты решения задачи для M =1,2 и 3 и точного решения для двух точек отрезка x =0 и x =1. Таблица 2.2 Значения температур для х=0 и х=1 на основе точного решения (ТР) и приближенных решений при М=1,2,3 X M=1 M=2 M=3 TP 0 1/3 -0,095 0,068 0 1 1/3 0,762 0,925 1 2.5. Естественные краевые условия Рассмотренный выше вариант МВН позволяет избавиться от необходимости подбирать для аппроксимации решений функции, удовлетворяющие краевым условиям. Однако при этом усложняется схема расчета, требуется больше параметров для достижения заданной точности решения, появляется необходимость вычисления интегралов вдоль границ, содержащих производные от искомой функции. Вычисление таких интегралов может оказаться затруднительным, если граничные поверхности криволинейны или имеют особенности. Для преодоления перечисленных недостатков можно преобразовать первое слагаемое в (2.24) к виду [1]: (2.26) где - линейные дифференциальные операторы более низкого порядка, чем D . При подстановке (2.26) в исходное уравнение (2.24) может оказаться, что часть второго слагаемого в (2.24) и последнее слагаемое в (2.26) отличаются лишь весовыми функциями. Поэтому, выбирая определенным образом весовые функции , можно добиться, чтобы эти слагаемые (содержащие интегралы от производных функций f по граничной поверхности Г ) взаимно уничтожались. Такая ситуация возможна лишь для определенных краевых условий, называемых естественными. Сама же формулировка задачи, основанная на преобразовании (2.26), носит название слабой формулировки МВН. В качестве примера рассмотрим решение задачи теплопроводности. при краевых условиях . Аппроксимацию искомой функции примем в виде (2.16) удовлетворяют условию T =0 при x =0, но не удовлетворяют условию . В качестве конкретных значений таких функций можно принять . Тогда в соответствии с выражением (2.24) уравнение МВН для рассматриваемой задачи примет вид . Первое слагаемое этого уравнения преобразуем согласно (2.26) Тогда уравнение МВН приобретет вид Выберем в качестве весовых функций Тогда уравнение МВН в окончательном виде может быть записано В этом уравнении не содержится слагаемых, содержащих производную на границе x =1, а краевое условие при x =1 является естественным. Параметры разрешающей системы алгебраических уравнений для рассматриваемой задачи принимают вид ; . Используя метод Галеркина Wn = Nn и выбирая базисные функции в виде Nm =xm можно получить: , . При М =2 из решения системы получаются следующие значения параметров am : a1 =11.758 a2 =3.458. Cоответствующие значения для точек x =0.5 и x =1.0, производной в точке x =1.0 и сравнения их с точными значениями приведены в таблице 2.3. Таблица 2.3 Результаты решения задачи на основе точного решения (ТР) и метода Галеркина (МГ) при М =2 Функция МГ ТР 6.743 6.754 15.216 15.232 18.67 20.0 2.6. Общая формулировка естественных краевых условий для задач теплопроводности Рассмотрим уравнение теплопроводности T =T (xi ), Q =Q (xi ) – заданная функция координат xi , при граничных условиях , ni – направляющие косинусы границы Г . Представим аппроксимацию решения этого уравнения в виде разложения по базисным функциям , (2.27) где . Для определения констант am используем общую формулировку уравнения МВН для =1,2… (2.28) Первое слагаемое преобразуем с помощью формулы Остроградского-Гаусса В результате уравнение МВН примет вид Примем также . При этом . Таким образом, краевое условие на Г q становится естественным, а уравнение МВН приобретает вид . (2.29) Подставляя в (2.29) представление в виде (2.27) получим [K ]{ }=[R ], , (2.30) . 2.7. Применение МВН в задачах теории упругости Рассмотрим применение МВН для решения задач теории упругости. Полный комплект уравнений, необходимых для решения задач теории упругости в перемещениях рассмотрен в первой главе и представлен соотношениями (1.1)−(1.21). В отличие от рассмотренных выше задач теплопроводности для скалярных функций , в задачах теории упругости искомая функция является векторной. Поэтому МВН в задачах теории упругости формулируется для системы уравнений, порядок которой определяется размерностью пространства решаемой задачи. Соответствующим образом скалярные базисные N m и весовые W n функции заменяются аналогичными векторными функциями и . . Рассмотренные в первой главе уравнения можно непосредственно использовать для решения задачи в рамках МВН, однако оказывается более удобным преобразовать их к виду, приводящему к слабой формулировке МВН. Запишем уравнения равновесия (1.12) и граниченые условия (1.1 и 1.13) для задач теории упругости в форме (2.24)МВН: для =1,2… (2.31) Первое слагаемое в этом уравнении можно привести к виду: Положим ; ; ; Тогда , т.е. условия на - естественные краевые условия. В результате уравнение МВН в слабой формулировке может быть записано n =1,2, . . . М (2.32) В матричной форме это уравнение можно представить в виде (2.33) ; В соответствии с обозначениями, принятыми в § 1.1 первой главы ; уравнение (2.33) примет вид (2.34) Конкретный вид матриц и приведен в § 1.1 первой главы. По аналогии с (2.27) положим ; (2.35) ; ; . Тогда уравнение МВН в окончательном виде может быть записано (2.36) Таким образом, уравнение МВН для задач теории упругости может быть сведено к системе алгебраических уравнений относительно параметров аппроксимации искомых функций : (2.37) ; , n, m =1,2... M где { }-вектор, составленный из параметров , конкретная структура которого, согласованная со структурой матрицы [K ] и вектора {R }, определяется из условия построения оптимального алгоритма решения системы (2.37). Представленная выше формулировка МВН может быть непосредственно получена на основе известного в теории упругости принципа возможной работы [2], согласно которому сумма работы, совершаемой внутренними и внешними силами в теле, находящемся в равновесии, при бесконечно малых виртуальных (не нарушающих кинематических граничных условий) перемещениях равно нулю: (2.38) В матричном виде уравнение (2.38) может быть записано в виде (2.39) Представим аппроксимации перемещений и в виде ; . В соответствии с обозначениями первой главы ; ; Тогда уравнение (2.39) примет вид для (2.40) В силу произвольности параметров из (2.40) следует уравнение (2.36), полученное на основе традиционной формулировки МВН . 3. Метод взвешенных невязок с кусочным определением базисных функций (метод конечных элементов) 3.1. Особенности задания базисных функций в методе конечных элементов (МКЭ) В предыдущих главах предполагалось, что при аппроксимации функций , базисные функции Nm определены на всей области V и интегралы (3.1) вычислялись сразу для всей области V . Рассмотрим другой возможный способ аппроксимации искомой функции f [1-3]. Для этого разделим область V на ряд непересекающихся подобластей (конечных элементов e =1-E ) и введем аппроксимацию функции f отдельно для каждой подобласти. Тогда при обеспечении необходимых условий гладкости таких функций вдоль границ конечных элементов (КЭ): ; (3.2) ; при условии, что , . Если подобласти имеют простую форму, а базисные функции определяются однотипно, то можно эффективно применять МВН для областей сложной формы. При этом все рассмотренные ранее процедуры МВН могут быть использованы для каждого отдельного конечного элемента.