Учебно-методическое пособие Рекомендовано методической комиссией механико-математического факультета для студентов ннгу, обучающихся по направлению подготовки 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

<

br />

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)

;

при условии, что , .

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

Таким образом, в МКЭ базисные функции задаются на локальных носителях – конечных элементах.

В МКЭ наряду с КЭ вводятся узлы, расположенные на границах, а иногда и внутри КЭ.

Все узлы КЭ области упорядочены, каждый узел характеризуется номером и связями с конечными элементами, содержащими такой узел.

Базисные функции в узлах задаются следующим образом:

Nm
=1 - в узле с номером m
,

Nm
=0 - в узлах n
m
и на КЭ, не содержащих узел m
. (3.3)

Nm
– отлична от нуля на КЭ, содержащих узел m
.

С учетом этих особенностей стандартная аппроксимация функции в МКЭ может быть записана в виде

, (3.4)

где , а – значение функции в узле m
. Иначе говоря, параметры приобретают в МКЭ конкретный физический смысл, а функция принимается равной нулю, так как граничное условие для f
может быть легко удовлетворено в граничных узлах

.

В представлении (3.4) базисные функции Nm
носят название “глобальных базисных функций”.

В свою очередь на каждом КЭ с индексом «е» глобальная аппроксимация (3.4) может быть выражена через значения функций в узлах i
-ого текущего КЭ и локальные базисные функции элемента (функции формы КЭ)

(3.5)

где n
– число узлов КЭ

Функции формы определены на соответствующих КЭ «e
» и принимают значения

в узле i
, в других узлах .

Обычно в качестве базисных функций в МКЭ используются степенные полиномы.

Например, аппроксимация функции одной переменной может быть представлена в виде полинома

, (3.6)

где – определяются из условия равенства функции в узлах узловым значениям функции .

Подставив полученные значения в уравнение (3.6) и сгруппировав слагаемые при можно получить аппроксимацию функции в виде (3.5):

.

Согласно такой интерпретации базисных функций в МКЭ аппроксимация (3.4) для линейных КЭ может быть представлена двояким образом.

Пусть требуется определить при

где и − глобальные базисные функции в узлах и , соответственно.

Эту же функцию можно представить через функцию формы КЭ

где и - локальные базисные функции элемента «r».

3.2. Аппроксимация решения дифференциальных уравнений с использованием кусочно-определенных базисных функций

Рассмотрим особенности применения кусочно-определенных базисных функций для аппроксимации решения дифференциальных уравнений в рамках МВН.

Пусть требуется определить некоторую функцию в области V
, удовлетворяющую следующей системе уравнений

,

.

В соответствии с общей схемой МВН аппроксимация решения поставленной задачи может быть осуществлена на основе решения уравнений

или

.

Рассмотрим условия, которым должны удовлетворять кусочно-определенные базисные функции приведенных выше уравнений в рамках МВН. Очевидно, эти условия должны обеспечить существование интегралов в исходных уравнениях, так как при кусочном задании базисных функций производные от таких функций могут претерпевать разрыв[1,3,4].

Согласно известным требованиям условия интегрируемости функции на сегменте сводятся к условию ее ограниченности на рассматриваемом сегменте и наличию ограниченного числа точек разрыва. Иначе говоря, кусочно-непрерывная и ограниченная на данном сегменте функция, интегрируемая на этом сегменте.

Применительно к рассматриваемой проблеме эти условия можно интерпретировать следующим образом: если операторы D
или G
содержат производные порядка “S
”, то базисные функции должны обеспечивать кусочную дифференцируемость функции до порядка S
-1.

Другими словами, базисные функции должны принадлежать классу гладкости : функция принадлежит классу , если она и все ее производные до порядка
r являются непрерывными.

Например, в случае интерполяции допустимы разрывные функции, так как S
=0.

Для рассмотренных выше уравнений теплопроводности и теории упругости S
=2, поэтому базисные функции должны принадлежать классу С
1
.

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

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

В качестве примера рассмотрим применение МКЭ для решения одномерной задачи теплопроводности.

,

при краевых условиях

.

Разделим отрезок на М
конечных элементов и построим аппроксимацию искомой функции в виде

, (3.7)

где - значение температуры в узле m
.

Уравнение МВН для этой задачи может быть представлено в виде

. ( 3.8)

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

В приведенном уравнении базисные функции должны принадлежать классу гладкости С
1
, то есть должны быть непрерывными и и первая производная , поэтому линейная аппроксимация для функций здесь не годится.

Для снижения требований гладкости можно преобразовать уравнение (3.8) к виду, соответствующему слабой формулировке МВН

(3.9)

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

, (3.10)

,

.

В свою очередь, при кусочном задании базовых функций

, (3.11)

где - вклад конечного элемента e в общую матрицу Knm
.

Нетрудно заметить, что ненулевой вклад будут давать только те элементы, которые содержат узлы с номерами n и m. При использовании линейной интерполяции отличными от нуля функциями на элементе «е
» будут Nm
и Nm
+1
,
причем .

Для определения вклада типового элемента введем локальную систему координат ,

где и - глобальные координаты первого и второго узлов элемента «e»,
h
– длина элемента.

В локальной системе линейные функции формы элемента (локальные базовые функции) могут быть записаны в виде

. (3.12)

При этом аппроксимация функции в элементе может быть представлена

,

где (i
=1−2)- значение в i-
ом узле элемента.

Для типового КЭ, содержащего узлы m
и n

,

или

.

Отсюда видно, что при вычислении нет необходимости производить вычисления для всех КЭ, достаточно произвести их для одного элемента в локальной системе координат. Если - функции формы КЭ, i
=(1,2) – номер локального узла элемента, то матрица [k
]e
будет иметь вид

Далее необходимо просуммировать вклад отдельных элементов в общую матрицу .

Рассмотрим процесс такого объединения для области, составленной из трех элементов (М
=3) одинаковой длины h
=1/3.

Для элемента e
=1 ненулевой вклад в knm
дают базовые функции с номерами n,
m
=1,2. В локальной системе этим номерам соответствуют

.

Для элементов с номерами e
=2,3 аналогичные вклады могут быть записаны в виде

,

.

Компоненты вектора {R
} будут равны

В результате система алгебраических уравнений (3.10) примет вид:

.

где (m
=1-4) –глобальные неизвестные в узлах сетки КЭ.

Вычеркивая из полученной системы строки 1 и 4 (поскольку и известны из краевых условий) и подставляя в строку 3 значение , получим

Решая систему, найдем .

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

Для сравнения в таблице 3.1 приведены результаты решения этой задачи, полученные на основе метода конечных разностей (МКР) и точного решения (ТР).

Таблица 3.1

Результаты решения задачи на основе точного решения (ТР), МКЭ и метода конечных разностей (МКР) при М
=3

Функция	МКЭ	МКР	ТР.
T (x =1/3)	0.2885	0.2893	0.2889
T (x =2/3)	0.6098	0.6107	0.6102
	0.8496		0.8509
	1.3156		1.3130

Сравнение полученных результатов позволяет судить о том, что для рассмотренной задачи при одинаковом числе разбиений отрезка точность решения на основе МКЭ и МКР практически одинакова.

В случае решения рассмотренной выше задачи с краевыми условиями при x
=0 и при x=1.в уравнение МВН должна быть включена невязка краевых условий на границе x
=1, где задан поток тепла

.

Выполняя интегрирование по частям и, полагая , получим

При использовании метода Галеркина и делении отрезка на три равных КЭ, матрицы [k
]e
отдельных элементов матрица [K
] всей системы сохраняются прежними. Меняется лишь вектор {R
}

С учетом решение задачи сводится к решению системы 3-х алгебраических уравнений. Результаты решения и сравнение их с результатами решений, полученных на основе точного решения (ТР) и метода конечных разностей (МКР) приведено в таблице 3.2.

Функция	МКЭ	МКР	ТР
	0.2193	0.2168	0.220
	0.4634	0.4576	0.4648
	0.7600	0.7493	0.7616
dT/dx (x =0)	0.6457		0.6481

Таблица 3.2

Результаты решения задачи на основе точного решения (ТР), МКЭ и метода конечных разностей (МКР) при М
=3

Представленные данные показывают, что на основе МКЭ результаты решения задачи оказываются значительно более точными, чем на основе МКР с использованием центрально-разностной аппроксимации градиента температуры на конце отрезка (погрешность ).

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

4. Построение базисных координатных функций в МКЭ

4.1. Основные требования к координатным функциям в МКЭ

Первым и вероятно наиболее ответственным шагом построения конкретных методик решения задач на основе МВН в форме МКЭ является выбор базисных функций используемых конечных элементов. Очевидно, что для получения качественного решения такие функции должны удовлетворять определенным требованиям, обеспечивающим сходимость КЭ решений [2-4]. Часть таких требований была сформулирована в предыдущей главе для базисных функций в МВН, а именно:

- функции должны быть линейно независимы;

-функции должны удовлетворять условиям полноты.

В МКЭ кусочно-определенные базисные функции также должны удовлетворять условиям непрерывности - иметь непрерывные производные до порядка S
-1, где S
- максимальный порядок производных, входящих в подынтегральные выражения МВН.

Кроме этого к координатным функциям в МКЭ предъявляется еще ряд специфических требований, среди которых важнейшими являются следующие:

- функции должны быть согласованы с формой КЭ, числом и расположением степеней свободы в элементе. Число степеней свободы элемента должно соответствовать числу неопределенных параметров в полиномиальном представлении функции;

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

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

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

4.2. Построение базисных функций конечных элементов в обобщенных координатах

Наиболее естественным образом базисные функции элементов могут быть получены на основе решения системы алгебраических уравнений для значений выбранного степенного полинома в узлах элемента. Этот подход является достаточно общим, применим для любых типов элементов и законов распределения функций и носит название «базисные функции в обобщенных координатах» [5]. Основные этапы применения этого подхода удобнее продемонстрировать на примере построения функций формы четырехузлового прямоугольного элемента (см. рис.4.1).

Закон изменения функции u
(x
,y
) в таком элементе может быть записан в виде степенного полинома

,

или

, (4.1)

{a
T
}=(a
0
a
1
a
2
a
3
)

Далее составляется система алгебраических уравнений для значений функции в узлах элемента

, (4.2)

,

Рис. 4.1. Схема четырехузлового прямоугольного элемента

Неопределенные параметры выражаются через узловые функции элемента

В результате функции формы элемента могут быть получены в виде

(4.3)

.

Одним из основных недостатков рассмотренного подхода является необходимость обращения матрицы для каждого конечного элемента исследуемой области. В случае использования полиномов выше первого порядка обращение этой матрицы в явном виде представляется весьма сложным и приходится применять численные процедуры. Кроме этого при использовании неполных полиномов выполнение условий межэлементной совместности функций, построенных на основе рассмотренного подхода, не является очевидным. Для выяснения этих условий следует записать закон распределения функций для каждой из сторон (граней) на основе принятого распределения в элементе и установить соответствие между числом неопределенных параметров в получающемся распределении с числом узловых значений функций на рассматриваемой стороне. Например, в представленном на рис.4.2 четырехугольном конечном элементе с билинейным законом распределения функций (4.1) непрерывность функции u
(x
,y
) вдоль сторон 1-2 и 2-3 будет выполняться по разному. Если в элементе ввести локальные координаты и направленные вдоль первой и второй сторон соответственно, то связь этих координат с общими декартовыми координатами x
,y
можно представить в виде

.

Рис. 4.2. Схема четырехузлового непрямоугольного элемента

Тогда изменения функции и вдоль сторон S1
и S2
могут быть определены подстановкой этих выражений в соотношения (4.1)

.

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

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

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

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

4.3. Построение локальных систем координат

Наиболее существенные недостатки рассмотренного выше общего метода построения базисных функций в МКЭ могут быть устранены путем получения этих функций в явном виде с использованием специально выбранных локальных систем координат [2,3,5]. Обычно эти координаты нормированы и связаны со сторонами элемента.

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

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

,

где - координаты центра элемента; -длина элемента.

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

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

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

В двумерном случае для прямоугольных элементов локальные координаты выбираются так, чтобы стороны прямоугольника совпадали с координатными линиями (см. рис. 4.3). Эти координаты связаны с глобальными координатами x
и y
соотношениями

(4.4)

где - координаты центра элемента, обозначенные на рис. 4.3; 2а и 2b – размеры сторон элемента вдоль осей x
и y
, соответственно.

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

Рис. 4.3. Локальная система координат четырехузлового прямоугольного элемента

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

(4.5)

где и - площади треугольников с вершинами в точке P
и основаниями и соответственно; - площадь элемента. Координаты и определяются аналогично. Основаниям треугольников соответствуют значения

Рис. 4.4. Локальная система координат треугольного элемента

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

,

, (4.6)

.

Решая (4.6) относительно , можно получить

, (4.7)

где i
, j
, k
образуют циклическую перестановку индексов 1.2.3.

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

4.4. Лагранжево семейство элементов

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

(4.8)

где - координаты узлов интерполяции.

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

Например, в локальной системе (- координата центра элемента, − его длина) функции формы элемента первого порядка будут иметь вид

,

(4.9)

Для квадратичного элемента функции формы будут такими

(4.10)

Функции формы прямоугольных двумерных элементов могут быть получены в виде произведения одномерных полиномов Лагранжа по каждой из координат и . При этом, в общем случае, порядок полиномов р
- вдоль и q
- вдоль координаты может быть различным. Если обозначить номера узлов в элементе индексами i j
, где i
и j
- номера узлов интерполяции на сторонах элемента (см. рис. 4.5) вдоль координат и соответственно, то сама искомая функция и функция формы элемента могут быть записаны в виде

. (4.11)

Рис. 4.5. Схема нумерации узлов элементов Лагранжева семейства

В частности, для билинейного элемента (см. рис. 4.3) p
= q
=2

(4.12)

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

. (4.13)

Рассмотренное семейство лагранжевых элементов удовлетворяет условию непрерывности функций вдоль граней и сторон элемента и условиям геометрической изотропии (при m
= n
). Основными недостатками элементов этого семейства являются наличие внутренних узлов в элементах, функции которых аппроксимируются полиномами выше первого порядка, и неполнота аппроксимирующего полинома. Например, характеристический полином, аппроксимирующий функцию в лагранжевом прямоугольном элементе 2-го порядка, может быть записан в виде

. (4.14)

Этот полином не содержит членов со степенями и , которые были бы желательны для полного кубического представления функции, и в то же время включает член более высокого порядка.

4.5. Сирендипово семейство элементов

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

Рис. 4.6. Характер изменения функций в сирендиповых Лагранжевых элементах

Как уже говорилось, функции формы сирендиповых прямоугольных элементов первого порядка совпадают с лагранжевыми и определяются соотношениями 4.12.

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

- угловые узлы:

; (4.15)

- на сторонах:

- на сторонах:

.

Функции формы элементов сирендипова семейства удовлетворяют условиям непрерывности вдоль граней и сторон элемента и условиям геометрической изотропии и также представляют функции в элементе в виде неполных полиномов порядка p
+1 (где p
- порядок полинома вдоль сторон элемента). Однако в отличие от лагранжевых элементов они не содержат членов более высокого порядка. В частности, характеристический полином для сирендипова элемента второго порядка будет содержать лишь первые восемь членов в выражении (4.14).

4.6. Треугольные элементы

Одним из существенных недостатков рассмотренных выше четырехугольных и шестигранных элементов является невозможность построения для них базисных функций, соответствующих представлению функции в виде полного полинома соответствующего порядка. Для треугольных и тетраэдральных элементов построение координатных функций в виде полных полиномов не вызывает каких-либо затруднений [2,3,5].

Действительно, рассмотрим семейство конечных элементов, полученных в результате деления сторон треугольника на m
равных частей при последовательном увеличении m
, начиная от m
=1. Каждому уровню такого разделения соответствует элемент порядка m
, имеющий узлы на пересечении сторон треугольника и линий, соединяющих точки деления сторон (см. рис. 4.7). Общее число таких узлов точно равно числу неопределенных параметров в полном полиноме степени m
.

Рис. 4.7. Схема иерархии треугольных конечных элементов

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

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

, (4.16)

,

(4.17)

Рис. 4.8. Схема нумерации узлов треугольных элементов различных порядков

Для элементов первого порядка (m
=1) функции формы совпадают с L
координатами

. (4.18)

Связь L
координат с глобальными координатами X
и Y
осуществляется на основе соотношения (4.7).

Для элементов второго порядка (m
=2) функции формы элемента принимают следующий вид:

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4.19)

.

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

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

и т.д.

При интегрировании по площади треугольника величин, зависящих от L
координат, может быть использовано следующее соотношение [5]

(4.20)

Все приведенные выше соотношения для треугольных элементов легко распространяются на пространственные элементы в виде тетраэдров.

4.7. Изопараметрические элементы и численное интегрирование

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

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

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

Идея использования функций формы для введения криволинейных координат принадлежит Тайгу, а ее развитие для широкого класса различных элементов связано с именами Айронса и Зенкевича. Отображение из локальной системы координат в декартову x,y,z
осуществляется посредством соотношений (см. рис.3.9)

(4.21)

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

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

Например, в локальной системе координат изменение функции в этом же элементе может быть записано в виде

(4.22)

Если порядок функции выше, чем порядок (соответственно m
>n
), то элемент носит название "суперпараметрический", если наоборот, то элемент - "субпараметрический". Чаще всего порядок функций и сами функции принимаются одинаковыми (m
=n
), и такие элементы носят название "изопараметрических".

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

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

В то же время непосредственная реализация изопараметрических элементов требует решения ряда дополнительных вопросов, которые не возникли при использовании рассмотренных ранее “простых элементов”.

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

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

В соответствии с (3.25) вычисление ее коэффициентов матрицы жесткости элемента сводится к интегрированию матрицы

. (4.23)

Элементарный объем в криволинейных координатах преобразуется по известной формуле

. (4.24)

Если криволинейные координаты являются нормализованными координатами, рассмотренными ранее в разделе 4.3, то интеграл (4.23) может быть записан в виде

. (4.25)

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

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

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

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

, (4.26)

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

(4.27)

При вычислении интегралов по двум и трем переменным, формулы интегрирования по схеме Гаусса имеют вид

, (4.28)

,

где - количество точек интегрирования по каждому из направлений.

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

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

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

P
= 1; n
= 2,

P
= 2; n
= 3,

P
= 3; n
= 4. (4.29)

Минимально допустимый порядок интегрирования определяется требованием, чтобы точно вычислялся объем конечного элемента [2,3], определяемый выражением (4.24). Использование более низкого порядка интегрирования может привести к нарушению сходимости решения по мере сгущения сетки элементов. В случае изопараметрических элементов, входящий в выражение (4.24) det
J
представляет собой полином от локальных координат , порядок которого легко установить, зная функции формы элемента. В частности, для рассмотренных элементов минимально допустимое число точек будет таким

P
= 1; n
= 1,

P
= 2; n
= 2, (4.30)

P
= 3; n
= 3.

При определении минимально допустимого порядка интегрирования следует также иметь в виду, что, согласно имеющимся исследованиям [6], максимальная скорость сходимости конечно-элементных решений сохраняется при точном интегрировании всех членов энергии деформации, соответствующих полному полиному наивысшего порядка, представленному в функции формы элемента. Минимально допустимое число точек согласно этому условию для рассмотренных выше элементов совпадает с (4.30).

Литература

1. Зенкевич О., Морган К. Конечные элементы и аппроксимация. М.: Мир, 1986. 318 с.

2. Капустин С.А. Метод конечных элементов в задачах механики деформируемых тел. Учебное пособие. Н.Новгород, 2002. 180 с.

3. Зенкевич О.К. Метод конечных элементов в технике. M.:Мир,1975. 544 с.

4. Стренг Г., Фикс Дж. Теория метода конечных элементов. M.:Мир, 1977. 349 с.

5. Норри Д., де Фриз Ж. Введение в метод конечных элементов. М.: Мир, 1981. 304с

6. Fried I. Numerical Integration in the Finite Element Method //Int. J.Comput. and Struct. 1974. V.4, N 5, p.921-932.