Решение задач тепло- и массопереноса с помощью метода точечных источников поля Текст научной статьи по специальности «Математика»

УДК 517.63:548.55

РЕШЕНИЕ ЗАДАЧ ТЕПЛО- И МАССОПЕРЕНОСА С ПОМОЩЬЮ МЕТОДА ТОЧЕЧНЫХ ИСТОЧНИКОВ ПОЛЯ

© 2006 г. С.Ю. Князев, Е.Е. Щербакова

Большое число задач тепло- и массопереноса сводится к численному решению дифференциальных уравнений в частных производных эллиптического и параболического типов. К таким уравнениям приводит также значительный ряд задач электромеханики, теории упругости, течения жидкостей и газов и т.д. Численное решение этих задач может быть затруднено, если поверхности (контуры), на которых задаются граничные условия, имеют сложную форму, особенно в том случае, когда граничные условия содержат нормальные производные от искомой функции. Использование метода точечных источников поля (МТИ) позволяет в ряде случаев значительно упростить нахождение численного решения задачи [1]. Наиболее результативен МТИ при решении краевых задач для однородного эллиптического уравнения с постоянными коэффициентами.

Однородное уравнение эллиптического типа

Известно [2], что всякое однородное уравнение эллиптического типа с постоянными коэффициентами может быть приведено к уравнению Гельмгольца

AU ± k 2U = 0.

(1)

Уравнение Лапласа будем рассматривать как частный случай уравнения Гельмгольца при k = 0. Фундаментальные решения трехмерного уравнения Гельмгольца, которые при знаке «плюс» в левой части (1) имеют вид

g (r ro ) =

exp (±ik |r - ro |)

r - r0

а при знаке «минус»

g (r o ) =

exp (±k |r - ro I)

r - r 0

(2)

(3)

g (r, ro ) = ln (1/|r - r o|).

В соответствии с принципом суперпозиции и физической постановкой решаемой задачи потенциал поля, описываемого уравнением (1), можно представить в виде

U (r ) = £ qg (r, R j ),

(4)

где радиус-вектор Я ^ определяет положение заряда q. Численное решение уравнения Гельмгольца с

помощью МТИ сводится к нахождению зарядов {д ^},

создающих поле, близкое к искомому.

Пусть требуется найти решение уравнения (1) в области V, ограниченной поверхностью (для двумерной задачи - контуром) Б, на которой выполняются условия

a(r )U (r )+ß(r)

dn

= Y(r)s .

(5)

Расположим заряды поля вокруг границы Б

на некотором расстоянии от нее (вне области V). Число зарядов N положим равным числу узловых точек на границе Б, радиус-векторы которых г. Подберем заряды д ^ таким образом, чтобы условия (5) выполнялись во всех узловых точках границы. Подставив (4) в (5) и записав полученные выражения для всех узловых точек границы Б, получим систему линейных уравнений для определения величин зарядов q ^.

Уравнение этой системы с номером i имеет вид

N

Е

j=i

(

a(ri )g(ri,R j )+ß(rг)

dg (r i, R j)

\

dn

q} =Y(ri). (6)

можно рассматривать как потенциал поля, созданного в точке г единичным положительным зарядом q = 1, помещенным в точку г0. При к = 0 формулы (2), (3) дают выражение для фундаментального решения уравнения Лапласа g (г, г0) = 1/ |г - г01. Для двумерного уравнения Гельмгольца фундаментальные решения выражаются через цилиндрические функции [1], а для двумерного уравнения Лапласа

После решения системы (6) и нахождения зарядов qполе в произвольной точке области V можно найти с помощью формулы (4).

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

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

Приведем в качестве иллюстрации возможностей МТИ численное решение двух тестовых задач. В первой задаче решались двумерные уравнения Лапласа и Гельмгольца для двухсвязной

V = {x2 + 4y 2 < 4}П j( x + +

(7)

Граничные условия имели вид (5) с постоянными значениями а и в, равными единице. В качестве точного решения Ц0(г) использовался потенциал поля, созданного несколькими точечными зарядами, величина и положение которых (вне области V) задавались с помощью генератора случайных чисел. Функция у(г) в правой части (5) соответствовала точному решению задачи. Критерием точности служили значения средней и максимальной относительной погрешности решения

(U(r)-U0 (r))

е = -

reV

тах

(UО (r ))| reV - min(Uо (r)|

(U(r)-Uо (r))|

тах (

(UО (r ))| reV - min (Uо (r))|

(8)

тах

Узлы на границе располагались равномерно; заряды д^ распределялись также равномерно вдоль

периметра границы на некотором расстоянии от нее. Часть зарядов помещалась внутри кругового отверстия, часть - за пределами эллипса. Вычисления показали, что значения средних и максимальных относительных погрешностей имели отличие менее одного порядка (обычно это отличие было не более 2-3 раз). Поэтому в дальнейшем в качестве меры погрешности используется средняя относительная погрешность е. На рис. 1 приведены результаты численных экспериментов.

зоо

200

Е

100

1о

1о

1о-

Ю"

1о-

о

-- 1о-

1о

Ю"6

ю-8

Из рис. 1 видно, что с ростом числа узлов п погрешность решения е быстро убывает до значений порядка 10-6 Такая точность достигается при относительно небольших значениях узлов (п ~ 200). Заметна зависимость точности решения от расстояния к, на котором заряды отстоят от границы 5. При небольших значениях к точность решения, как и следовало ожидать, достаточно низкая (кривая 1 на рис. 1, полученная при к = 0,05). При увеличении к точность решения возрастает (кривая 2 на рис. 1, полученная при к = 0,1; кривая 3 - при к = 0,3). Максимальная точность достигается при к, сравнимых с линейными размерами области V. При дальнейшем росте к точность практически не меняется. При значительных к (в несколько раз превышающих линейные размеры области V) наблюдается небольшой рост погрешности вычислений.

Во второй задаче решалось трехмерное уравнение Гельмгольца ДЦ - к 2Ц = 0 для области, ограниченной сферой радиуса Я. Граничные условия (5) брались при постоянных значениях параметров а = в = у = 1. При таких условиях уравнение Гельмгольца имеет аналитическое решение. С ним сравнивался результат численного решения, и по формулам (8) оценивалась относительная погрешность. На рис. 1 приведены зависимости относительной погрешности Е от числа узлов Ы, качественно согласующиеся с результатами первой тестовой задачи. При оптимальном расположении зарядов (на расстоянии, равном к = 2Я от поверхности сферы) максимальная точность достигается уже при числе узлов N = 500 и имеет порядок 10-15 (кривая 4 на рис. 1). Кривые 5 и 6 на рис. 1 получены при к = 0,5Я и к = 0,2Я соответственно.

Неоднородные уравнения эллиптического типа

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

AU ± к 2U = f (r).

(9)

Рис. 1. Зависимость относительной погрешности от числа узлов при решении уравнения Лапласа для двухсвязной области (кривые 1-3) и Гельмгольца для сферы (кривые 4 - 6)

Уравнение Пуассона будем рассматривать как частный случай неоднородного уравнения Гельмгольца при к = 0. Пусть для уравнения (9) известно также частное решение Ц0 (г) (не обязательно с однородными условиями на границе):

ДЦ 0 (г )± к 2Ц 0 (г ) = / (г).

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

/(г) = ¡1 (*)+ /2 (У), или в виде / (г) = ¡з (г). В°з-можны также другие варианты, когда можно получить частное решение уравнения Гельмгольца в аналитическом виде. Во всяком случае, всегда можно найти близкую к /(г) функцию ф(г), для которой удается

получить частное решение и0 (г) уравнения (9) и оценить допускаемую при такой замене ошибку.

и

n

S

Введем функцию u (r) = U (r) - U0 (r). Эта функция удовлетворяет однородному уравнению Гельмгольца Au ± к 2u = 0 с граничными условиями

а—+ ßu dn

dU 0

= 1 Y-а-

dn

-ßU о

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

Уравнения параболического типа

Точечными источниками можно моделировать не только стационарные поля, удовлетворяющие эллиптическим уравнениям Гельмгольца и Лапласа, но и нестационарное поле, описываемое уравнением параболического типа, например уравнением диффузии или теплопроводности. Фундаментальное решение уравнения теплопроводности [1]

ф(, r0, t, t о ) =

1

f

(a (t -1 о ))

d/2

exp

r - r 0

\

4D (t -10 )

(10)

таким образом, чтобы граничные условия (5) выполнялись во всех узловых точках во все моменты времени I ^. Заряды qij находим для каждого временного

уровня последовательно, начиная с уровня . Предположим, что все заряды для уровня 1п -1 найдены. Введем обозначения:

N n—1

1 (r ) = С 0 + ЕЕ Яг; Ф ( R i, ПТ, Л) >

i =1 j =1

Яг = Ягп .

(12)

Для нахождения зарядов qi = qin на временном уровне / п подставим (11), с учетом (12) в граничные условия (5). В результате получаем систему уравнений для определения величин зарядов qi. Уравнение этой системы с номером т имеет вид:

N

Е|_а mn ф(Г m , R г > nT (n - !)Т) +

i =1

+ ß mn "дф- (rm > Rг , nZ, (n - 1)т)

Яг =

'(rm )-ßm

du (rm )

dn

можно рассматривать как потенциал поля, созданного в точке r в момент времени t единичным мгновенным источником, возникшим в момент времени 10 в точке с координатой r0. a - коэффициент температуропроводности; величина d равна размерности задачи (d = 2 или d = 3 ). Рассмотрим нахождение численного решения уравнения диффузии в объеме V при начальном условии: U(r,0) = C0 = const. В граничных условиях (5) теперь допускается зависимость параметров а, в и у не только от координаты точки на границе S, но и от времени. Разместим на границе S узловые точки. Введем шаг по времени т. В точках R i, расположенных на некотором расстоянии от границы S и число которых N равно числу узловых точек на S, поместим серию мгновенных зарядов q^,

последовательно возникающих в моменты времени tj = jT, где j = 1,2,3.... Поле, создаваемое системой

зарядов qtj в точке r, в момент времени t имеет

потенциал

N n

U (r, t ) = с 0 + EE q4 <p(r, R i, t, jT), (11)

i=1 j=1

где n < t / т .

Так как все слагаемые в правой части выражения (11) являются решениями уравнения диффузии, то отсюда следует, что и U (r, t) также удовлетворяет уравнению диффузии. Остается подобрать заряды qi]-

где а тп =а(гт , 1п ) Р тп = Р(гт , 1п ) У тп = У(гт , 1п )

В качестве примера приведем численное решение двух тестовых задач. Первой является задача для двухсвязной области (7), решением которой является точечный мгновенный источник (10). Положение мгновенного источника внутри кругового отверстия области (7) задавалось с помощью генератора случайных чисел. Граничные условия имели вид (5) с а = в = 1. у(г) в правой части (5) соответствовала точному решению задачи. Второй решалась задача об остывании цилиндра, имеющая точное аналитическое решение, выражаемое через функции Бесселя [2]. Варьировались: число узлов N (равное числу точечных мгновенных источников поля), расстояние к, на котором располагались точечные источники от границы Б области (7) и шаг по времени т. Для двух рассматриваемых задач обнаружилось близкое соответствие зависимостей относительной погрешности численного решения от указанных параметров. Устойчивость численного решения наблюдалась, если шаг по времени находился в интервале 2-10-3 < т/т0 < 20-10-3, где т0 = Я2/а -характерное время для данной задачи, Я - длина малой полуоси эллипса для первой или радиус цилиндра для второй задачи. Если шаг по времени находится в указанном интервале, то, как видно из рис. 2 погрешность решения сначала убывает со временем, а затем начинает возрастать. Причем, если шаг по времени оказывался вблизи границы устойчивости, то рост ошибки со временем был весьма быстрым (кривая 1 на рис. 2).

10-

10-

10-

3 4

2

1

0

0,05

0,1

0,15

thn

Рис. 2. Зависимость относительной погрешности от времени: т/т0 = 310-3, 5-10-3, 10-2 и 10-2 для кривых 1, 2, 3 и 4 соответственно; h = 0,1R (кривые 1-3), h = 0,3R (кривая 4)

Аналогичная зависимость наблюдается для расстояния h, на котором располагаются точечные заряды. Устойчивость численного решения наблюдалась, если расстояние h находилось в интервале 0,005 < h/R < 0,35. Если расстояние h оказывалось вблизи границы устойчивого решения, то относительная ошибка, начиная с некоторого момента, быстро возрастала со временем (кривая 4 на рис. 2). Заметная зависимость погрешности от числа узлов наблюдается только при их небольшом количестве. При увеличении числа узлов от 20 до 100 погрешность уменьшается более чем на 2 порядка от 3-10-2 до 2-10-4. При дальнейшем росте числа узлов погрешность практически не изменяется (рис. 3, кривая 1). Это, видимо, связано с тем, что при N > 100 численная погрешность определяется конечностью шага по времени т, которая при увеличении числа узлов не меняется.

t/t0 120 100

10

10

10

-80

60

40

20 0

0

50

100

150 N, Nt

Рис. 3. Зависимость относительной погрешности 8 от числа узлов N (кривая 1) и зависимость времени счета ^ от числа шагов по времени N1 (кривая 2)

Согласно формуле (11) с ростом времени и, соответственно, с увеличением числа шагов по времени число слагаемых в ряде (11) увеличивается и, следовательно, увеличивается объем машинного времени, затрачиваемого на вычисление искомого потенциала и (г, /). Кривая 2 на рис. 3 отображает зависимость машинного времени (в относительных единицах),

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

Совместное использование МКР и МТИ

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

Рассмотрим для определенности решение уравнения теплопроводности для области V, ограниченной поверхностью X Допускается зависимость коэффициента теплопроводности от температуры и, следовательно, нелинейность уравнения теплопроводности. Введем в области V, разностную сетку с шагом к. Узлы разностной сетки, у которых расстояние до границы вдоль какой-либо линии сетки не превышает к, назовем приграничными узлами. Прочие узлы, находящиеся внутри области V, назовем внутренними узлами. Граничными узлами будем считать точки пересечения линий разностной сетки с граничной поверхностью 5. Проведем через приграничные узлы замкнутую поверхность 5. Заключенную между поверхностями 5 и 5 часть области V обозначим как 5 V. Очевидно, что время установления теплового равновесия, время релаксации в подобласти 5V намного меньше времени релаксации во всем объеме V и это время сокращается с уменьшением шага разностной сетки к. Если шаг к достаточно мал, то время релаксации в подобласти 5V может оказаться много меньше используемого шага по времени т. В этом случае, пренебрегая изменением коэффициента температуропроводности по ширине зазора между 5 и 5, температурное поле в подобласти 5V приближенно можно описать уравнением Лапласа. Следовательно, на каждом шаге по времени значение температурного поля во внутренних узлах можно находить с помощью МКР, а в подобласти 5V - с помощью МТИ. Для этого

у

у

внутри области V на некотором расстоянии от поверхности 5 и вне области V на некотором расстоянии от граничной поверхности Б располагаем точечные заряды, число которых равно суммарному числу граничных и приграничных узлов. Температурное поле в подобласти 8V представляем в виде ряда (4). Решение задачи будет состоять в получении системы линейных уравнений для значений температуры во внутренних узлах разностной сетки и точечных зарядов внутри и вне области V. Каждому узлу разностной сетки (внутреннему, приграничному и граничному) будет соответствовать одно линейное алгебраическое уравнение. Для внутренних и приграничных узлов эти уравнения получают с помощью разностной аппроксимации уравнения теплопроводности. При этом значения поля в приграничных и граничных узлах, входящие в разностные уравнения, представляются в виде ряда (4). Для граничных узлов с помощью граничных условий получаются уравнения, аналогичные (6). После решения полученной СЛАУ с помощью ряда (4) вычисляют температуры в приграничных и граничных узлах, затем переходят к вычислениям для следующего шага по времени.

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

О (т) + О (к2). Установлено, что смешанный метод

позволяет получить решение задачи с удовлетворительной точностью. На рис. 4 приведены зависимости средней относительной погрешности от шага разностной сетки, полученные при двух различных значения шага по времени, равных т = 0,01-к2/а (кривая 1) и т2 = 0,1-к2/а (кривая 2); а - коэффициент температуро-

проводности. Для кривой 1 эта зависимость аппроксимируется выражением е ~ 0,3-й1'6.

s _

410-2

2E0-2 -

0,1 0,2 0,3 h/R

Рис. 4. . Зависимость относительной погрешности от шага по координате h

Таким образом, точность смешанного метода оказалась близкой к точности используемой разностной схемы.

Литература

1. Алексидзе М.А. Фуедаментальные функции в приближенных решениях граничных завдач. - М.: Наука, 1991. -352 с.

2. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М., 1977.

3. Баттерфилд Р., Бенерджи П. Методы граничных элементов в прикладных науках. М., 1984.

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

5. Бахвалов Ю.А., Гречихин В.В., Юфанова Ю.В. Комбинированные модели в расчетах электромагнитных полей // Изв. АН. Серия физическая. 2004. Т. 68. № 7. С. 1019-1022.

Южно-Российский государственный технический университет

(Новочеркасский политехнический институт) 19 июня 2006 г.

УДК 666. 71. 2

ВЛИЯНИЕ ФАЗОВОГО СОСТАВА И СТРУКТУРЫ ДЕКОРАТИВНОЙ СТЕКЛОМОЗАИЧНОЙ ПЛИТКИ НА ЕЁ СВОЙСТВА

© 2006 г. Е.А. Лазарева, А.М. Напрасник, Л.В. Дьяченко, В.В. Кирюшенко

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

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