Posobie_MathCAD_v2

		i x j	x		b	b	j	x i x dx
	aij
					a	a					.
									b
	i b j b i a j a									j	x i x dx
									a

Как видно из определения базисных функций, на концах от-
резка a, b	все они обращаются в ноль. Поэтому два первых
члена в этом выражении равны нулю и
b	n	xk 1
aij	j x i x dx	j x i x dx .
a	k 1 x		k

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

тах). Т.е. j x 0 при i 1 j i 1 . Отсюда следует вывод, что полученная нами матрица A является трехдиагональной. Элементы главной диагонали (i=j) определяются как

n

xk 1

2

xi

1

xi 1

1

2

k

1

i dx

h

2

dx

h

aii

dx

x

k

x

i 1

x

i

.

xi

xi 1

1

2

dx dx

2

h

h

xi

xi 1

На соседних с главной диагональю

n

xk 1

xi 1

1

1

dxdx

1

xi 1

1

a

dx

dx

,

i,i 1

2

i

i

1

h

h

k 1

x

k

x

h

h

x

i

i

n

xk 1

xi

1

1

1

a

dx

dx

.

i,i 1

i

i

1

h h

h

k 1

xk

xi 1

Таким образом, матрицу A можно записать как

131

	2		h	1	h		0			0		...			0
	1	h		2		h	1	h		0		...			0
A	0			1	h		2		h	...		0			0			.
	0			0			...			1	h	2		h	1	h
	0			0			...			0		1	h		2		h
Матрица коэффициентов этой													системы называется также
матрицей жесткости.							Систему								можно решить точным
												Ac

экономичным методом прогонки, описанным в главе 2. В других случаях возможно применение итерационных методов решения СЛАУ.

В случае если краевые условия имеют более сложный вид, необходимо ввести замену переменных с тем, чтобы привести их к виду u(a)=u(b)=0. Более подробно о методе, его вариациях, обобщении на многомерные задачи можно узнать из литературы

[4, 6].

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

1.Разбить область на подобласти (в общем случае, это могут быть неравномерные трехмерные элементы);

2.Выбрать вид базисной функции;

3.Для каждого элемента построить функционал и применить процедуру минимизации. Таким образом, определяются локальные коэффициенты СЛАУ;

4.Объединить локальные матрицы каждого элемента в глобальную матрицу жесткости (ансамблирование) и создать вектор правых частей;

5.Учесть граничные условия;

6.Решить трехдиагональную СЛАУ численно.

Пример 6.6. Рассмотрим задачу из примера 6.5, но примем, что наша нагрузка линейно изменяется в пространстве. Предположим, что Q есть линейная функция u. Для простоты положим Q=u. В качестве граничных условий возьмем u(0)=0 и u(1)=1.

132

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

1.Разбиваем отрезок [0,1] на M равных частей (формируем массив x)

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

иN2, соответственно)

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

4.В этой же ПФ происходит ансамблирование коэффициентов

исоздание матрицы жесткости K, которая является выход-

ным параметром работы ansmb. Далее создается вектор правых частей f.

5.В процедуру-функцию solv пересылаются количество узлов, вектор правых частей, граничные условия

6.Здесь же решается трехдиагональная СЛАУ методом прогонки (см. главу 2) и результат записывается в массив.

Невязка полученного решения для M=3 приведена на рис. 6.17. Можно убедиться, что порядок ошибки мал по сравнению с решением.

133

M 3

i 1 M 1

a 0

b 1

h

b

a

0.333

0 (i 1) h

M

xi

fi1 0

fiM 1

N1(i y)

xi 1 y xi

N2(i y)

y

xi 1

xi 1

xi

xi

ansmb(M)

for

i 1 M

for j 1 M

Ki j 0

for

i 1 M 1

hj xi 1

xi

hj

d

d

ke

2

N1(i y)

N2(i y) N1(i y) N2(i y) dy

1

dy

dy

0

ke2 1 ke1 2

hj

2

d

2

dy

ke

N1(i y)

(N1(i y))

1

1

dy

0

ke2 2 ke1 1

for

m 1 2

for

l 1 2

Ki 1 m i 1 l Ki 1 m i 1 l kem l

K

KK ansmb(M 1)

Рис. 6.17. Листинг решения примера 6.6. методом конечных элементов (начало)

134

f 1

1

f

M 1

1

xM 1

xM

x2

x1

i 2 M

f i 0

solv(M KK FF fl fr )

L M

for

i 2 L 1

CCi KKi i

BBi KKi i 1

AA i KKi i 1

3

KK2 3

KK2 2

3

FF2

KK2 2

for

i 3 L 1

BBi

i 1 CCi i AA i

AA i i FFi

i 1 CCi i AA i fi L fr

for i L 1 1

fi i i 1 fi i 1 i 1

fi

a1 solv(M 1 KK f fi1 fiM)

Рис. 6.17. Листинг решения примера 6.6 методом конечных элементов (окончание)

135

Глава 7. Решение уравнений в частных производных

7.1. Основные понятия. Постановка задачи

Для описания реального физического процесса (прогиба балки под действием нагрузки, движения газа в некотором объеме, распространение электромагнитных волн и пр.) строятся физи- ко-математические модели, на основе которых можно анализировать процессы качественно и количественно. Задачи описания движения сплошных сред (газа, жидкости, твердых тел), а также задачи теплопроводности, теории упругости, электрических и магнитных полей, и многие другие приводят к дифференциальным уравнениям. Независимыми переменными в физических задачах являются время t и пространственные координаты x, y, z. В качестве зависимых переменных в разных задачах используются компоненты скорости частиц среды, плотность, давление, температура, упругие напряжения, деформации и др. характеристики.

Допустим, что требуется найти решение на временном промежутке [t0, t1] в некоторой области изменения независимых переменных G(x,y,z). Математическая постановка задачи состоит из дифференциального уравнения (уравнений), а также дополнительных условий, позволяющих выделить единственное решение среди семейства всех решений данного уравнения. Дополнительные условия, заданные при t=t0 называются начальными данными, а условия, заданные на границе области G(x,y,z) – граничными или краевыми условиями. В качестве начальных и краевых условий, как правило, задают значения искомых функций и их производных. Задачи, у которых имеются только начальные условия, называются задачей Коши. Задачи с начальными данными и граничными условиями называют смешанной краевой задачей или нестационарной краевой задачей.

При исследовании установившихся состояний или стационарных (не зависящих от времени) процессов используются уравнения, не зависящие от времени. В этом случае решение ищется в области G(x,y,z), на границе которой задаются граничные условия. Такие задачи называются краевыми.

136

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

7.1.1. Характеристики. Типы уравнений.

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

Autt Butx Cuxx Dut Eux F .

Здесь u(t, x) – искомая функция, t, x – независимые переменные, A, B, C, D, E и F – коэффициенты уравнения, которые, вообще говоря, могут зависеть от t, x и u.

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

137

Если коэффициенты A, B, С – нулевые, а D 0 и E 0, то уравнение имеет первый порядок и называется уравнением переноса (адвекции). Если хотя бы один из коэффициентов A, B, С отличен от нуля, уравнение имеет второй порядок и может быть классифицировано по типам аналогично кривым второго порядка. Классификация уравнений связана с наличием характеристик – особых направлений, вдоль которых исходное уравнение может быть записано в виде полного дифференциала, и, следовательно, может быть проинтегрировано. Уравнение характеристик для уравнения в частных производных второго порядка имеет вид

dx		B	B 2 4 AC		.
dt			2 A

Количество характеристик зависит от знака дискриминанта B2 – 4 A C. Если он положителен, то уравнение (7.1) имеет две вещественных характеристики и называется гиперболическим. В случае нулевого дискриминанта уравнение имеет одну вещественную характеристику и является параболическим. Эллиптические уравнения не имеют вещественных характеристик (дискриминант отрицателен).

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

7.1.2. Аппроксимация и устойчивость

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

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

138

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

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

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

Lu , u0 v0 , lu ,

с помощью конечно-разностной схемы:

Lhuh h , u0 h v0 , lhuh h .

Рассмотрим также конечно-разностную задачу с «возмущенными» правыми частями, начальными данными и краевыми ус-

ловиями:
~ ~ ~0	~	~ ~
Lhuh h , u	v0	, lhuh h .

Говорят, что решение разностной задачи непрерывно зависит от входных данных задачи, если существуют такие константы M1>0, M2>0, M3>0, не зависящие от t и h, что

139

~

M1

~

M 2

~

M 3

~

(7.0)

uh (t) uh (t)

h h

h h

v0 v0

Для линейных разностных схем доказано, что устойчивость по начальным данным эквивалентна устойчивости по правой части, т.е. для линейных операторов Lh, lh достаточно показать, что

~

(t)

M

v

~

.

uh (t) uh

0 v0

Если независимых переменных несколько (например, x и t), то вводится понятие условной и безусловной устойчивости. Устойчивость называется безусловной, если соотношение (7.0) имеет место при произвольном соотношении между шагами разностной сетки и h. Схема будет условно устойчивой, если для выполнения (7.0) шаги по независимым переменным должны подчиняться дополнительным соотношениям.

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

Подробнее о свойствах аппроксимации и устойчивости будет рассказано ниже при рассмотрении примеров уравнений в частных производных.

7.2. Параболические уравнения 7.2.1. Постановка задачи

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

	u	2u	(7.1)
	t A	x2 F (x.t).

Здесь А>0 – константа (коэффициент теплопроводности), u(x,t) – искомое решение, F(x,t) – правая часть. Будем искать

140