Конечно-элементное моделирование локальных газодинамических процессов в трехмерных пористых структурах Текст научной статьи по специальности «Математика»

УДК 539.31

Ю. И. Д и м и т р и е н к о, А. И. Л е в и н а, А. А. Г а л и ц ы н

КОНЕЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ ЛОКАЛЬНЫХ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ В ТРЕХМЕРНЫХ ПОРИСТЫХ СТРУКТУРАХ

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

E-mail: dimit.bmstu@mail.com

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

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

Асимптотический анализ уравнений газовой динамики для пористой среды. Следуя работам [3-7], рассмотрим пористую среду периодической структуры, поры которой заполнены линейно-вязким совершенным газом (рис.1). Обозначим l0 — линейный размер ячейки периодичности У^ среды и x0 — характерный размер всей области V среды. Полагая, что выполняется соотношение: е = l0/x0 << 1, где е — малый параметр, введем локальные и глобальные безразмерные коор-

Рис. 1. Модель трехмерной пористой периодической структуры

динаты: ^ = 1x4 г и = х7х0. Тогда все функции, описывающие течение газа в порах, можно считать квазипериодическими; дифференцирование этих функций осуществляется с помощью следующего

д/ 1 д/ д/ „ правила: -— = — т—- + -т—-. Далее используются следующие обозна-дх, 8 дд, дх,

чения: / = —: / = —, / = —, V / = —. Область, занятую од-

7/г д£. '' д.Х/ Л дГ ^ дх,.

ной порой, будем обозначать V* , а ее границу . Операция осредне-

/ \ J

J щ.

ния функций по ячейке периодичности имеет вид (/) =

где ( =| dV^g — пористость среды.

Изотермическое движение газа в пористой среде, занимающей область V, в рамках сделанных выше допущений описывается системой уравнений Навье - Стокса, которая в декартовых координатах имеет следующий вид:

Р + V, (pv,) = 0, , eV,

j-((pv,), (pvv1)) = -Pt v^., ,e V

vi y

Eu Re

(i)

= 0, p = Ap; i, j = 1,2,3,

здесь р, р, V. — безразмерные давление, плотность и вектор скорости газа, р0, р0, v0 — характерные значения функций; ^ — безразмерное время; А = Яр0в0/р0 — безразмерная константа, в0 — температура, Я — газовая постоянная; Ей = р0/р^02 — число Эйлера; Яе = р^0х0/1 — число Рейнольдса, где ¡1 — вязкость газа.

Введем соотношения между тремя безразмерными параметрами задачи е, Ей, Яе следующим образом:

Eu0=Eusk ,Re0=Resm

(2)

где к и т — некоторые целые числа, которые соответствуют различным режимам течения газа в порах. Значения Ей0, Яе0 для определенности полагаем близкими к единице: Ей0 = 0(1), Яе0 = 0(1). Подставляем соотношения (2) в систему (1), получаем систему с одним малым параметром е:

Р +

t+s h )/г+(pvi), г=0 (pvi )t+s

s

pv.v . i J

V J J

+

1 s

= -S p/i-p i+-

k + m

/J

f

f

pv.v .

iJ

V J J

J

Eu0Re0

11

— v. , .. + — v. , .. .2 i / JJ s i / JJ

(3)

.. + v.

i,JJ

v

i

Z

0, p = Ap.

^sg

Осредняя функции p, p, v получим

(pvt) = p, (p) = p (p) = p

(4)

где р, р, V. — средние значения функций. Соотношения (4) в зависимости от режима течения газа либо служат дополнительными условиями, предъявляемыми к локальным параметрам течения в У^ либо являются обозначениями для вычисленных осредненных функций.

Согласно асимптотическому методу осреднения [7], решение системы (3) находим в виде асимптотических разложений:

v, = v, (0)( x ,£■, t) + sv,11>( x , t) + О (s2),

0)/

p = p(0)( x ,£■, t) + sp{i,( xt , t) + О (sO

o)/

(5)

p = p(0)( x ,£., t) + sp (1)( x л,, t) + 0(s2).

o)/

Подставим разложения (5) в систему (4) и сгруппируем члены при одинаковых степенях е. В зависимости от значений к и т возможны семь различных случаев:

к = 0, т = 1; к = 0, т = 2;

к = 1, т = 0; к = 1, т = 1; к = 1, т = 2. (6) к = 2, т = 0; к = 2, т = 1;

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

Формулировка локальной задачи. Далее будем рассматривать случай к = 2, т = 0. Случаи к = 0, т = 1 и к = 0, т = 2 были рассмотрены в [7, 8].

Для этого случая давление р(0)(х, ,$) = р(х, )и плотность р(0)(х ,$) =

= р (х ) не зависят от локальных координат В, Приравнивая нулю члены при в уравнении неразрывности и при е° в уравнении движения системы (3), получаем локальную задачу «на ячейке периодичности» (ЯП) для вычисления функций V р(1), р(1):

V,(0) = 0, $ еУея, -Р/,(1) + Ц% = Р,Л $ *Уея V,(0) = 0, $ , (7)

(р(1)) = 0,

где ц0 = —0—0. системе (7) добавим условия периодичности фун-

Еи Ке

кций, обозначаемые [[/]] = 0, здесь -1/2 < В < 1/2 — область ячейки периодичности У^ . Градиент давления р(г0) в системе (7) зависит только от глобальных координат х , и рассматривается как «входные данные» локальной задачи переноса (7).

Будем рассматривать далее трехмерную пористую структуру, у которой поры образуют продольные каналы вдоль всех трех осей ОВ (рис. 1). В силу линейности локальной задачи (7) ее решение можно представить в виде линейной функции от входных данных, т. е. от р,(0):

р(1)=£р(а)($)р(о), V(1)=Л £щ)р,<°>, (8)

а=1 Ц а=1

где функции Р(а)(£), Щг(а)(Вг) зависят от В

Подставив выражения (8) в локальную задачу (7), после исключения градиента р(г0) получим набор локальных задач для определения функций Р(а)(В ), Ща)(В ), которые в отличие от задачи (7) не содержат

констант, описывающих физические свойства газов, и не зависят от входных данных:

wa = 0,

-P™ +AWi(а) = h(a), 6 eV4 w(а) = 0, 6 GZ,sg,

(9) (10) (11) (12) (13)

(P(a)) = 0, [[W(a)]] = 0, [[P(a)]] = 0.

Здесь АЩ(а) = Щ^. Решение задач (9)-(13) определяется только внутренней геометрией пор, поэтому их решение применимо для расчетов фильтрации любых газов и жидкостей в рамках сделанных выше допущений.

В уравнениях (10) обозначены функции

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

Метод решения локальной задачи для ЭБ-пористых структур.

Для 3Б-структур в задаче (9)-(13) индекс г пробегает значения г = 1, 2, 3. Для каждого фиксированного значения а задача (9)-(13) представляет собой стационарную задачу течения некоторой фиктивной линейно-вязкой несжимаемой жидкости. Эта задача имеет ряд специфических свойств. Во-первых, поскольку функция р(1) в выражении (5) представляет собой пульсацию давления р по отношению к среднему значению р (из разложений (5) следует ер(1) = р - р(1) + 0(е2), а (р(0)) = р в силу (р(1)) = 0), то функция р(1) может быть как положительной, так и отрицательной. Во-вторых, наличие условия (12) делает задачу (9)-(13) интегрально-дифференциальной, что вместе с условиями периодичности (13) принципиально отличает ее от классической задачи Стокса.

Упростить решение задачи (9)-(13) позволяет следующая теорема о продолжении решения, являющаяся аналогом теоремы из работы [8].

Теорема. Пусть ячейка периодичности У^ ЗБ-структуры имеет зеркальную симметрию относительно координатных плоскостей и тогда решение Р(а, задачи (9)-(13) полу-

чим с помощью симметричного или антисимметричного продолже-

10 при i Ф а или пр и i = а и Pао) = 0,

И{а) = \ ,а

(14)

1 пр и i = а и Ра(а0) Ф 0,

Рис. 2. Перемена знаков у функций при симметричном и антисимметричном продолжении решения задачи (15)

ния (согласно правилу, показанному на рис. 2) функций Р(а), Щ.(а), определенных в 1/8 ячейки периодичности (в первом квадранте

У$§ = {$ : 0 < $ < У}) и являющихся решениями следующих задач:

Щ7(? ) = 0,

-Р(«) + ЩИ = $ еУ$„, (15)

Ща = 0, $ е^.

На рис. 2 показан способ симметричного и антисимметричного продолжения функций Р(1), Щ(1), во всю ЯП. Для остальных функций Р(а), Щ.а), а = 2, 3, решение продолжается аналогичным образом, путем циклической перестановки индексов и функций. Знак «+» на рис. 2 означает, что при переходе из первого квадранта функция не меняет знак (симметричное продолжение), а знак «—» означает, что меняет (антисимметричное продолжение). Легко убедиться, что конфигурация симметричных и антисимметричных продолжений функций р'(1), удовлетворяет полной системе (9)—(13). Доказательство теоремы осуществляется тем же методом, что и в работе [9].

Используя симметричное и антисимметричное продолжение решения системы (15), можно записать граничные условия на граничных плоскостях 1/8 ЯП У :

£ = 0 8 £ = 1: ^ = 0; Щ^2(1) = 0; Щ" = 0; Р(1) = 0;

2

1 дй(1) дй(1) дР(1)

62 = 0 8 62 = 1: = 0; й™ = 0; = 0; дР- = 0; (16)

£ = 0 8 £ = 1: ^ = 0; ^ = 0, ^ = 0, = 0.

3 3 2 д63

Вариационная формулировка локальной задачи. Введем обозначения для компонент тензора напряжений Т(а) и тензора скоростей деформаций В(а):

Т(а) = -Р (а)5г + Б(а), Ба) = ^(Щг + 1). (17)

У У У У 2 ' 1 11

Тогда уравнение равновесия в задаче (15) можно представить в дивергентном виде: Та = . Если ввести кинематически возможное поле скоростей 8 Ща как векторное поле, удовлетворяющее граничным условиям 1-го рода на части поверхности ЯП Ща) = 0, £ е ,

и вариацию давления 8Р(а), то, умножая скалярно уравнение равновесия в (15) на 8Щ(а), а уравнение несжимаемости на 8Р(а) и интегрируя получившиеся выражения по области У получим систему двух вариационных уравнений:

| Ща)8Р{а)ёУ8 = 0,

Уж (18)

| Т^5йгаёУг = | 4а)8Щ, £ еУ*.

У У

ж ж

Воспользуемся формулой преобразования произведения векторов из [10]:

Т^ЗЩа) = ТаЩ(а) )/1 - Т!ата. (18)

Тогда, используя формулу Остроградского - Гаусса, систему (18) можно записать следующим образом:

Г Б^дР(а)ёУж = 0,

J г г ж

Уж (19)

Г Т^дБ^У = Г 8ГШг(а)^-| (а)<Уя, £ еУ*.

Уж 2ж Уж

Здесь использованы формулы (17) и введено обозначение для компонент вектора усилий Sг(а) = п^а на поверхности 2ж области Уж с ,

где п. компоненты вектора нормали к поверхности Заметим, что на всей границе области течения Vв соответствии с граничными условиями из задачи (15) работа вектора поверхностных усилий обращается в нуль. Система (19) образует вариационную постановку задачи (15).

Метод МКЭ для решения вариационной задачи. Для решения системы вариационных уравнений (19) был применен метод конечных элементов с 10-узловым конечным элементом в форме тетраэдра, который отличается от классического 10-узлового КЭ [11,12] тем, что имеет 34 степени свободы: по три компоненты скорости Жга) в каждом узле и по одному значению давления Р(а) в каждой вершине тетраэдра (рис. 3). Аппроксимация функций в каждом КЭ по скоростям была квадратичной, а по давлению — линейной:

{W} = [фШ

3 3x30 30

Р = {фр f {У},

1 4 4

(20)

где {Ж} = (Ж1(а), Ж2(а), Ж3(а)) — координатный столбец скоростей в КЭ;

— координатный столбец скоростей в узлах; [у] — координатный столбец давлений р = Р(а) в вершинах КЭ.

Матрица функций формы [Ф] и столбец {Ф } имеют следующий

вид:

[Ф] =

3x30

ф 0 0 ф 2 0 0 . ф ■■ ^10 0 0

0 ф 0 0 ф 2 0 . .. 0 ф 0

0 0 ф\ 0 0 ф 2 .. 0 0 ф ^10

где

{Фр г = (А, ь2, L3, L4),

ф = (2L - 1)Ц,Ф2 = (2L2 - = (2L3 - 1)4Ф = (2L4 -1)L4,

Ф5 = 4LiL2,06 = 4ЦЦ,фп = 4^4, ф, = 4LiL4,09 = 4^,ф0 = 4L3L4.

L = ¿T ( + brX + СУ + dZ ), 1 = 12,3,4 6V

6V =

1 X1 У1 Z1

1 X2 У 2 Z2

1 X3 У3 Z3

1 X4 У 4 Z4

(21)

a a , b , с , d получаются циклической перестановкой индексов 1, 2, 3, 4.

Координатный столбец скоростей деформаций [Б}Т = (Б11, Б22, Б33, Б12, Б13, Б23) представим в виде [Б} = [I]{Ж} = [!][Ф][^} = [5][?}, где матрица [Ь] имеет вид

д 0 0 1 д 2 д^2 1 д 2 д^з 0

0 д 0 1 д 0 1 д

2 д£ 2 д£

0 0 д д£з 0 1 д 2 д^з 1 д 2 д^2

(22)

Вариации компонент вектора скоростей и тензора скоростей деформаций, а также давленияр с помощью (8) представляются в виде:

S{W} = [Ф]5[д],

Sp = {ФрГ S{y}, S{D} = [B]S{q}.

(22)

Введем координатный столбец напряжений [Т}т = (Т^, Т^, Т^, Т1(2), Т(3), Т2(3а)), который на основании определяющих соотношений (17) связан с координатным столбцом скоростей деформаций линейным соотношением:

[Т} = -[Фр ][ у}+[Б}, (23)

6 6X4 4 6 V )

где [Фр] — матрица, образованная диадным умножением столбцов:

[Фр ] = [/с}[Фр }Т, а [/о}Т = (1,1,1,0,0,0), [1Х } = (1,0,0).

p 6x4

Введем также координатный столбец внешних усилий [£}Т = (£1(а), ^2(а), ^3(а)) на поверхности тела. Тогда, подставляя (20), (22), (23) в систему (19), получим систему линейных алгебраических уравнений (СЛАУ) для каждого отдельного КЭ:

[ K ]{q}- [ Bp ]{y} = {f }, [Bv ]T {q} = 0,

(24)

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

[K ] = J((B ]T [B])dVg, [Bp ] = J([B ]T [Фр ])),

vg Vg

Ш = J(( tf})) {fv } = J(C {h})dVg, f = Л + fv,

(25)

g

g

где {h} = (ha),h3a)). При формировании глобальной СЛАУ суммарные векторы правык частей {f.}, соответствующие вектору поверхностный усилий {5}, обращаются в нуль, поскольку все внутренние стороны КЭ при суммировании перечисляются дважды, но знак у вектора {5} при этом меняется из-за изменения знака вектора нормали на этих сторонах. На тех же сторонах КЭ, которые принадлежат внешней границе области V^ , работа вектора {5} обращается в нуль в силу граничных условий задачи, поэтому в глобальной СЛАУ {f,} = 0.

Поскольку матрица итоговой СЛАУ (24) несимметричная и сильно разреженная, возник вопрос о выборе численного метода ее решения, так как большинство алгоритмов, разработанных для решения разреженный матриц, предназначено для симметричнык матриц. В этой связи для решения СЛАУ были выбраны методы бисопряженных градиентов BiCG и TFQMR (Transpose-Free Quasi-Minimal Residue) [13]. Метод бисопряженных градиентов BiCG имеет по сравнению с методами FOM и GMRES гораздо более простую вычислительную схему и меньшее требование к памяти, однако BiCG, как и любой другой метод, использующий операции с транспонированной матрицей A, плохо поддается реализации на многопроцессорный вычислительных системах с распределенной памятью.

Был разработан программный комплекс, реализующий решение локальной задачи (19) описанным выше методом МКЭ. Комплекс решает следующие основные подзадачи: формирование расчетной области; задание свойств материалов; формирование на основе расчетной области сетки КЭ; решение задачи МКЭ; обработка результатов применения процедуры МКЭ.

Формирование геометрии расчетной области производилось в пакете SolidWorks. Для генерации конечно-элементной сетки применялся

Рис. 3. 10-узловой КЭ с тремя компонентами скорости в каждом узле и одним значением давления в вершинных узлах КЭ

генератор ИЭ. Визуализация расчетов проводилась при помощи программы Ме8ИУ1е^^, разработки кафедры ФН-11 МГТУ им. Н.Э. Баумана.

Тестовая задача: сравнение с аналитическим решением задачи для течения Пуазейля. Рассмотрим пористую одноканальную структуру, у которой имеются сквозные цилиндрические каналы в направлении ^ (см. рис. 4). В этом случае задача (15) имеет точное решение согласно [4]:

#1(1) =-(а2 -г2)/2, (26)

что соответствует течению Пуазейля для линейно-вязкой жидкости, где г2 = ^22 + £32; а — радиус канала. Расчет производился для случая а = 0,4.

Рис. 4. Сравнение численного МКЭ-реш ения задачи Пуазейля с аналитическим. Распределение безразмерной скорости в канале

На рис. 4 показаны результаты численного МКЭ-решения задачи (15) изложенным выше методом и аналитическое решение. Видно, что численный метод обеспечивает удовлетворительную точность вычислений. Численные данные сравнения приведены также в табл. 1.

Таблица 1

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

г 0 0,04 0,08 0,12 0,16 0,2

W -102 МКЭ 7,93 7,8 7,53 7,18 6,65 5,76

W -102 аналит 8 7,92 7,68 7,28 6,72 6

Г 0,24 0,28 0,32 0,36 0,4

W -102 МКЭ 4,95 3,84 2,59 1,11 0

W -102 аналит 5,12 4,08 2,88 1,52 0

Результаты расчетов для 3В пористых структур. Было проведено численное моделирование трехмерных пористых структур с одним сквозным каналом (рис. 5) и с тремя сквозными цилиндрическими каналами (рис. 6), которые соединяются в центре координат ЯП со сферической порой. В расчетах варьировалась геометрия пор, при этом исследовалось влияние относительного радиуса каналов на характер течения газа в порах.

Исследованные варианты геометрии приведены в табл. 2. В этой же таблице приведены максимальные значения скорости Ж1(1) и давления Р(1) газа в порах. Максимальное значение скорости Ж1(1) достигается в цилиндрической части пор на координатной оси О^, а максимальное значение давления — на твердой стенке на границе стыка сферической и цилиндрической частей пор (см. рис. 5). С увеличением радиуса как цилиндрической, так и сферической части пор максимальное значение скорости Ж1(1) возрастает, а давление Р(1) с ростом цилиндрической части сначала растет, а затем уменьшается.

Поперечные компоненты скорости Ж2(1) и Ж3(1) имеют локальные максимумы значений на границе стыка сферической и цилиндрической частей пор. Абсолютные значения этих компонент для рассмотренных геометрий примерно на порядок меньше максимальных значений скорости Ж1(1).

Таблица 2

Исследованные варианты геометрии 3D пористой структуры

Радиус сферической части поры R Радиус цилиндрических каналов г Максимальная скорость Г 04 О2 Максимальное давление /"11

0,1 0,05 0,135 0,052

0,2 0,05 0,177 0,115

0,3 0,05 0,249 0,180

0,3 0,1 0,830 0,187

0,3 0,15 1,575 0,233

0,3 0,2 2,242 0,197

0,3 0,25 3,358 0,1652

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

Рис. 5 (начало). Картины полей скоростей (а—в) и давления (г) в пористой однока-нальной структуре с Я = 0,3 и г = 0,1

а

б

в

Рис. 5 (окончание)

б

Рис. 6 (начало). Картины полей скоростей (а—в) и давления (г) в пористой структуре с тремя каналами, Я = 0,3 и г = 0,1

г

а

Рис. 6 (окончание)

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

в

г

СПИСОК ЛИТЕРАТУРЫ

1. С а н ч е с-П а л е н с и я Э.Теория колебаний и неоднородные среды. - М.: Мир, 1984. - 472 с.

2. Б а х в а л о в Н. С., П а н а с е н к о Г. П. Осреднение процессов в периодических средах. - М.: Наука, 1984. - 352 с.

3. Д и м и т р и е н к о Ю. И. Механика композиционных материалов при высоких температурах. - М.: Машиностроение, 1997. - 356 с.

4. Д и м и т р и е н к о Ю. И., Г л а з и к о в М. Л. Численный расчет проницаемости и процессов фильтрации в пористых средах. В сб. «Аэрокосмические технологии». - М.: Изд-во МГТУ им. Н.Э. Баумана. - 2002. С. 132-137.

5. Д и м и т р и е н к о Ю. И., Г л а з и к о в М. Л. Моделирование процессов фильтрации в периодических пористых средах. - Вестник МГТУ им. Н.Э. Баумана. Естественные науки. - № 1. - 2003. - С. 59-71.

6. Д и м и т р и е н к о Ю. И., Л е в и н а А. И., Б о ж е н и к П. Конечно-элементное моделирование локальных процессов переноса в пористых средах // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2008. - № 3. -С. 90-104.

7. Д и м и т р и е н к о Ю. И., И в а н о в М. Ю. Моделирование нелинейных динамических процессов переноса в пористых средах // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - № 1. - 2008. - С. 39-56.

8. D i m i t r i e n k o Y u. I. Dynamic Transport Phenomena in Porous Polymer Materials Under Impulse Thermal Effects// Journal of Transport in Porous Media. - Vol. 35. -N 2. - 1999.

9. Д и м и т р и е н к о Ю. И., К а ш к а р о в А. И. Расчет эффективных характеристик композитов с периодической структурой методом конечных элементов // Вестник МГТУ им. Н.Э. Баумана. Естественные науки. - № 2. - 2002. - С. 95-108.

10. Д и м и т р и е н к о Ю. И. Тензорное исчисление. - М.: Высшая школа. - 2001. -576 с.

11. С е г е р л и н д Л. Применение метода конечных элементов. - М.: Мир, 1979.

12. Ш а й д у р о в В. В. Многосеточные методы конечных элементов. - М.: Наука, 1989.

13. Y o u s e f S a a d. Iterative methods for sparse linear systems. - 2000.

Статья поступила в редакцию 27.10.2011.

Димитриенко Юрий Иванович родился в 1962 г., окончил в 1984 г. МГУ им. М.В. Ломоносова. Д-р физ.-мат.наук, профессор, заведующий кафедрой «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана, действительный член академии инженерных наук. Автор более 200 научных работ в области механики сплошной среды, вычислительной механики, термомеханики композитов, математического моделирования в науке о материалах.

Левина Александра Игоревна. Окончила МВТУ им. Н.Э. Баумана в 2003 г. Старший преподаватель кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор 8 научных работ. Область научных интересов: теория фильтрации, моделирование процессов в периодических структурах. Галицын Алексей Андреевич. Окончил кафедру «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана в 2011 г. В настоящее время — инженер ОАО «ВПК «НПО машиностроения».Область научных интересов: метод конечного элемента, методы оптимизации, теория фильтрации.