MATLAB.Exponenta
MATLAB и Simulink на русском
Технологии разработки и отладки
		сложных технических систем

Математика\Partial Differential Equations Toolbox

В.Е.Шмелев "Partial Differential Equations Toolbox. Инструментарий решения дифференциальных уравнений в частных производных"
1. Конечноэлементная технология решения задач математической физики

  В оглавление книги \ К следующему разделу \ К предыдущему разделу

1.2. Проекционная формулировка конечноэлементных методов

Существует вариационная и проекционная формулировка конечноэлементных методов. PDE Toolbox Matlab поддерживает только проекционную формулировку. Проекционные конечноэлементные методы являются частным случаем методов взвешенных невязок.

1.2.1. Математическая основа проекционных методов

Пусть в N-мерной области G евклидова пространства EN необходимо рассчитать распределение некоторой физической величины, которое описывается дифференциальным уравнением в частных производных (PDE):

L uf = 0 , x im1.gif (174 bytes) G                                            (1.3)

с граничными условиями

LГ u = f Г , x im1.gif (174 bytes) Г ,

где u=u(x) – анализируемая физическая величина; x – точка пространства EN ,
L и LГ – дифференциальные операторы; Г – граница области G; функции f, f Г описывают распределение внутренних и граничных источников исследуемого поля.

Рассмотрим пространство M   функций, определённых на области G, такое, что если u1(xim1.gif (174 bytes) M , u2(x)   im1.gif (174 bytes) M , то
u3(x) = a1u1(x)+a2u2(x)   im1.gif (174 bytes) M (любая линейная комбинация функций, принадлежащих пространству M, также принадлежит пространству M). Введём скалярное произведение функций:

(u,v) = dG ,

тогда нормой функции u I M будем называть неотрицательное число
||u|| = (u,u)1/2 , расстоянием между функциями u и v будем называть величину, равную норме ||uv||. Такое пространство M называется гильбертовым пространством [5].

Методы взвешенных невязок (проекционные методы) [1, 6] можно разделить на три большие группы:
1) граничные методы (если приближённое решение точно удовлетворяет дифференциальному уравнению);
2) внутренние методы (когда приближённое решение точно удовлетворяет граничным условиям);
3) смешанные методы (не удовлетворяются ни граничные условия, ни дифференциальное уравнение).
К граничным методам можно отнести метод Трефца, методы граничных элементов и др. [5, 7, 8, 9].

Рассмотрим некоторые внутренние методы. Введём невязку приближённого решения ua дифференциального уравнения (1.3): R = L uaf . Приближённое решение ua представим в виде:

ua(x) = u0(x) + ,

где – известные аналитические функции, которые часто называют пробными; ai – коэффициенты, которые подлежат определению; u0(x) выбирается так, чтобы функция ua по возможности наиболее точно удовлетворяла граничным условиям. Уравнения для определения коэффициентов ai можно получить исходя из условия:
(R,wk(x))=0, где k = 1, …, N’ ;
wk(x) называется весовой, или поверочной функцией. В соответствии с выбором вида функций wk(x) различают метод подобластей, метод коллокаций, метод наименьших квадратов, метод Галёркина (когда в качестве весовых функций берутся пробные: wk(x)=) и другие [6].

1.2.2. Метод Галёркина с конечными элементами

Важное место среди вариационных и проекционных методов занимают конечноэлементные методы (методы конечных элементов (МКЭ)) [1, 2, 3, 6]. Рассмотрим краевую задачу (1.3) для трёхмерной односвязной области V с границей S. Пусть область V+S пространственно дискретизована: будем считать, что она “натянута” на M узлов с фиксированными координатами xi, yi, zi, i=1, … M. Тогда конечноэлементная аппроксимация приближённого решения ua краевой задачи (1.3) имеет вид (1.1).

Пусть поверхность S состоит из двух частей: S=S1+S2. Пусть на части граничной поверхности S1 задано распределение искомой величины u:

u = поверхностное распределение,                                   (1.4)

а на части граничной поверхности S2 задано распределение нормальной или тангенциальной составляющей некоторого дифференциального оператора распределения искомой величины u:

n*LГ u или n**LГ u = поверхностное распределение.           (1.5)

(1.4) называются граничными условиями первого рода (условиями Дирихле), (1.5) – граничными условиями второго рода (условиями Неймана). При решении PDE (1.3) граничные условия первого рода удовлетворяются путём задания соответствующих значений величины u в граничных узлах.

Если использовать метод Галёркина применительно к аппроксимации (1.1), то дифференциальное уравнение (1.3) преобразуется в интегральное матричное уравнение относительно узловых значений искомой величины:

L ([N]*[u(у)]) dV = f dV .                           (1.6)

Матрица [u(у)] может быть вынесена за знак объёмного интеграла справа. Если дифференциальный оператор L содержит порядка z , то необходимо гарантировать кусочную дифференцируемость производных порядка z –1 от используемых функций формы [N] [1], иначе в точках их разрыва производные порядка z будут неограничены, что приведёт к появлению неопределённостей в интегралах, содержащихся в (1.6).

В [1] показано, что интеграл в левой части (1.6) может быть заменён выражением вида

*(D [N]) dV*[u(у)] + dS ,

где операторы E и D имеют более низкий порядок дифференцирования, чем исходный оператор L. Такая замена даёт возможность понизить требования к гладкости функций формы [N].

Модифицированное таким образом уравнение (1.6) можно преобразовать в алгебраическое матричное уравнение вида [K]*[u(у)] = [F] . Для выполнения такого преобразования интегралы в модифицированном уравнении (1.6) можно вычислять путём суммирования по конечным элементам. Однако поверхностные интегралы по внутренним границам элементов в сумме должны обращаться в нуль (если в расчётной области отсутствуют поверхностные источники моделируемого поля), иначе данное интегральное преобразование не будет иметь смысла. Если функции формы имеют полиномиальный вид, то объёмные интегралы могут быть вычислены аналитически [10 и др.]. Однако для элементов высокого порядка, как правило, используется численное интегрирование. Аналитическое интегрирование приводит к существенной экономии вычислительных ресурсов, однако на практике мало применяется из-за трудности получения расчётных формул и усложнения программирования.

Как и в конечно-разностных методах, основным свойством получаемых систем алгебраических уравнений является ленточная или профильная структура матрицы коэффициентов (матрицы “жёсткости” [2]), а также слабая заполненность внутри ленты или профиля. Поэтому для хранения таких систем в памяти ЭВМ используют технологию разреженных матриц [11, 12, 13]. Прямым методам технологии разреженных матриц свойственно заполнение матриц в процессе факторизации. Поэтому наряду с прямыми применяются также итерационные методы. Прямые и итерационные методы реализованы в ядре MATLAB в виде операций и функций над разреженными матрицами.

1.2.3. Применение метода Галёркина с конечными элементами к эллиптическому PDE

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

–div(c*grad(u)) + a*u = f .                                                     (1.7)

Для простоты будем полагать, что (1.7) записано относительно скалярного поля u, коэффициенты c и a также скалярные. Пусть коэффициенты c и a распределены в расчётной области неравномерно. Применение метода Галёркина с конечными элементами к уравнению (1.7) даёт следующее интегральное уравнение:

(–div(c*grad(u)) + a*u) dV = f dV ,             (1.8)

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

div([N]T*c*grad(u)) dV + grad[N]T*c*grad(u) + [N]T a*u) dV = f dV .

Применим теорему о дивергенции к первому интегралу последнего уравнения, перенесём его в правую часть и применим конечноэлементную аппроксимацию (1.1) к искомому полю u:

grad[N]T*c*grad[N] + [N]T a [N]) dV*[u(у)] = f dV + *c*grad(u)*dS

или короче [K]*[u(у)] = [F] ,                                   (1.9)
где [K] = grad[N]T*c*grad[N] + [N]T a [N]) dV – глобальная матрица коэффициентов (матрица “жёсткости” [2]);
[F] = f dV + *c*grad(u)*dS – глобальная матрица-столбец правой части (матрица “нагрузок”).

При выполнении расчётов глобальное матричное уравнение (1.9) формируется путём накопления при рассмотрении всех конечных элементов:
[K] = ; [F] = ,
где [K(e)] =grad[N(e)]T*c*grad[N(e)] + [N(e)]T a [N(e)]) dV – локальная (элементная) матрица коэффициентов;
[F(e)] = f dV – элементная матрица правой части;
[Fe)] =*c*grad(u)*dS – граничноэлементная матрица правой части;
V(e) – объём конечного элемента;
S(e) – поверхностный элемент границы расчётной области,
dS – векторный дифференциал площади граничной поверхности (понятие о векторе площади поверхности см. в [14]). Если распределение скалярных коэффициентов c, a и источника поля f  можно аппроксимировать кусочно-однородным (постоянным в пределах каждого элемента) распределением, то c, a, f   можно выносить за знаки соответствующих элементных интегралов.

Поверхностный интеграл в правой части (1.9) соответствует граничным условиям Неймана (1.5)
n
*c*grad(u) = поверхностное распределение.
При реальных расчётах эти граничные условия могут быть учтены путём вычисления данного поверхностного интеграла при рассмотрении соответствующих граничных элементов (граничными элементами называются грани конечных элементов, принадлежащие границе расчётной области). Граничные условия Дирихле (1.4) учитываются путём задания соответствующих значений величины u в граничных узлах. Задание значений в некоторых узлах конечноэлементной сетки осуществляется путём модификации матриц [K] и [F] [2]. При расчётах такая модификация производится в последнюю очередь непосредственно перед решением матричного уравнения (1.9).

При моделировании плоскопараллельных полей объёмные интегралы в (1.9) можно заменить на двумерные, а поверхностные – на криволинейные, учитывая, что
dV = z dU, dS = z dl**1z ,
и z можно сократить в левой и правой части (1.9). Здесь dU – дифференциал двумерной расчётной области в плоскости z=const, dl – векторный дифференциал длины границы расчётной области. Как было сказано выше, PDE Toolbox MATLAB поддерживает только симплекс–элементы и только двумерные краевые задачи. В связи с этим формулы для вычисления элементных матриц, соответствующих эллиптическому PDE (1.7), имеют вид:
[K(e)] = 0,5 c(e) (e)*([ax]T [ax] + [ay]T [ay]) + a(e) ;
[F(e)] = f(e) ;
[Fe)] = = *|le)|* ,       (1.10)
где c(e), a(e), f(e) – постоянные значения параметров уравнения (1.7) в пределах конечного элемента;
(e) – определитель элемента, равный его удвоенной площади;
le) = xe)*1x + ye)*1y – вектор длины граничного элемента;
be) = *1x + *1y = c*grad(u) – постоянное значение вектора в пределах граничного элемента (в соответствии с заданными граничными условиями Неймана).

Формулы (1.10) следуют из соотношений, свойственных двумерным симплекс-элементам [2]:
; ,
где L – отрезок прямой, l – его длина; a, b, c – произвольные целые числа >=0;
(e) = det .

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

  В оглавление книги \ К следующему разделу \ К предыдущему разделу


Поиск по сайту:

Система Orphus

Яндекс.Метрика