КОНСЕРВАТИВНА КIНЦЕВО-РIЗНИЦЕВА СХЕМА ЗАДАЧI СТЕФАНА ДЛЯ РIВНЯННЯ ДИФУЗIї Текст научной статьи по специальности «Математика»

УДК 519.6

КОНСЕРВАТИВНА К1НЦЕВО-Р13НИЦЕВА СХЕМА ЗАДАЧ1 СТЕФАНА ДЛЯ Р1ВНЯННЯ ДИФУЗИ © Славко Г.В.

кременчуцький державний пол1техн1чний унюерситет

им. Михайла Остроградського

факультет електрон1ки i комп'ютерно! 1нженерп

вул. Першотравнева, 20, м. Кременчук, 39600, Украша e-mail: emathiSmail.ru

Abstract. Conservative finit-difference solution's scheme of an equation of diffusion with a free boundary is adapted for a software realization in MathCAD. Developed program allows to research solutions for arbitrary function which determines a free boundary.

Метод кшцевих р1зниць - найбшын розповсюджений числовий метод розв'язання рпшянь з чаетинними похщпими, Клаеичним прикладом е заетоеування цього методу для розв'язання р1вняння дифузп - одновим1рного по проеторовш координат! рпшяння парабо. гршого типу. Задач1 криетал1зацп внмагають побудовн алгоритму розв'язання таких рпшянь. за умови вшьно'1 меж1, що визначаетьея деякою функ-щею, та граничних умов третього роду [1]. Специфжа мов програмування еучаеннх математичних програм Мар1е, Ма1Ь( 'а<1. Ма11аЬ та пшшх потребуе адаптацп вщомих а.норн ! ,\Гп; к1нцево-р1зницевих схем у вщповщноета з Тх семантикою з урахуванням оеобливоетей накоппчення похибки обчпелень [2]. Безперечш переваги заетоеування цих програм визначаютьея швидгаетю програмування та наочшетю подання отри-маних результате 1 компенеують вптратн часу на розробку таких алгоритм1в. Тим не менше безпоеередне програмування у цих програмах в бшыноета випадгав вияв-ляетьея бшьш ефектпвнпм, шж внкорпетання вбудованпх функщй.

Оеоблнв1етю задач1 (1)-(4) е вшьна межа, що визначаетьея функщею та гра-ннчн*1 умови третього роду (комбшащя шукано! функцп та И похщно! по проеторов1й координат!) за наявноета похщпо1 йз/сИ функцп, що визначае вшьну граннцю. Таким чином, маемо задачу типу Стефана:

Вступ

1. Задача Стефана у загальному виглядг

(1.1)

граничш умови:

(1.2)

гуЩЬИ + Su (L, t) = ipL(t),x = L,t> О,

(1.3)

початков! умови:

и (х, 0) = ф (х), s (0) < X < L, t = 0.

(1.4)

Граничш умови ¡ндн'п дозволяють побудувати явну кшцево-р1зницеву схему, 1Т перевагою безумовно е можливють отримання значень шукано! егаковоТ функцп на наетупному часовому шар1 через вщом1 значения попереднього шару, а тому прог-рамна реал1зацп проста (без розв'язання систем лшшнпх р1внянь), Але, як вщомо,

застосування ще! схеми обмежуеться вимогою умовоТ стшкоста о2ту/ (/гк)2 < 1/2 , яка

визиачае можлиш кроков1 характеристики си?кп т, /гк, 11аяншп ь вшьно1 меж1 обумо-влюе змшу для часових крогав просторових крогав, а це в свою черту накладае бшын жорп к*1 вимоги до сп кн. Таким чином, явну р1зницеву схему за критер1ями «нжноп *1 <-. п. г вважати неприйнятною для розв'язання ще! задач*!. Тому будуватимемо неявну схему, яка хоч 1 призводить до необхщноета розв'язання снстеми л'нппннх рпшянь. але е абсолютно етшкою, що особливо важливо для змшно! си?кп, Зазиачимо тут, що для неявно! схеми чисельний розв'язок у раз*! його зростання е завищеним, а у ра и спадания - занижении, Це вщбуваетьея за рахунок того, що числовий розв'язок за неявною схемою внзначаеться значениями егаково! функцп на верхньому часовому шарь

2, побудова консервативно! р13ницев01 схеми

Для побудови р1зиицево1 схеми для задач1 (1)-(4) скористаемося методикою робота [3], Побудуемо для облаета

* (/) < .г < I.. 0 < / < 1.1 < х:

кшшт,о-р1 шпиону сп ку

гк = к-т, /,• = оЛ7..у = с

М/гт = {%) = 8к ■ Ък, //' = (/- - 8к)/М, 8к = 8 } . Наведемо деяга пояснения. Крок за часом т е незмшним для ипх шар'п; с'п-ки т = I' / М. де М - к*1.1ьк*1с1 ь часових п перин, пи. а просторовий крок ¡гк змшюеться у ш. пюш. шоп *1 з виразом ¡гк = (1 — -чк)/АГ, де N - гальгасть просторових пперна. пи. '!.м'ша просторового кроку на кожному часовому крощ призводить до залежноета просторових координат хк ну ¡. пи си?кп вщ часового шару. Часов1 шари визначають-ся виразом I1' = к ■ т. Таким чином, будемо мати сп ку. яка адаитуеться до ¡.\iiini меЖ1 часово-просторово! область У вузлах сп?кп визначатимемо шукану функщю задач1 ик = и . Диференцшш операторп апрокспмуемо вщношенням кшцевнх

р1зннць для к -го I а к — 1 часових шар1в:

к+1 к+1 _ к

~ 3 3 +0(г), (2.1)

ди

т

д2

и

дх2

к+1 к+1 _ , к+1

= 3+1 )2+ + О ({¡гк+1у ) . (2.2)

Шдетавимо (5),(6) в (1), отримаемо для внутриншх вузл1в сп?кп:

ук+1 - ик ик+1 - 2ик+1 4- ик+1 / оч

3 - 3+1 3 + / (*к+1) + О (г + (^+1)2) . (2.3)

т

(Ьк+гУ

Якщо псшдш першого порядку, що ирису! н*1 в граничних умовах апрокеимувати вщношеннями правих 1. пипх кшцевпх р1зниць:

ди дх

к+1

3=0

и

к+1

и:

к+1

¡1к+1

+ О (кк+1)

ди дх

к+1

и

к+1 N

и

к+1 N-1

¡1к+1

+ О (кк+1) , (2.4)

то апрокеимащя буде наближепням першого порядку 1 загальний порядок апрокеи-мацп схемп буде першим, не зважаючи на те, що в шших вузлах порядок апрокепма-цп по проеторовнм змшннм другпй. Таким чином, для збереження другого порядку апрокеимацп, в граничних вузлах розкладемо и'\ * 1 в око. 1*1 точки х = з (I1' 1) в ряд Тейлора за змшною х до третьоТ похщно']" включно, а и'р"^ в такий же ряд в окол1 точки х = Ь , За умови, що функщя и (х. I} в граничних вузлах мае перни пох*1.111*1 за часом 1 друп по х , матпмемо:

ик+1

и(вк+1 + кк+1,гк+1)

и

к+1 N-1

и

и

и

к+1 N

к+1

3 (9) 1 (10) отримаемо:

ди дх

ди дх

к+1

0 к+1

N

и

к+1

и

+

ди дх

к+1

ди дх

к+1

N

к+1

¡1к+1 +

0

¡1к+1 +

д2

и,

дх2

д2и к+1 ^к+1

дх2 2

к+1

N 2 '

(2.5)

Iгк+1

и

к+1 N

и

к+1 N-1

1з (1) матимемо:

кк+1

д2и дх2

+

д2и к+1 кк+1

дх2 0 2

д2и к+1 Ьк+1

дх2 N 2

+ 0

+ 0

1 (ди „, Д

Переходячи в (13) до граничних вузл1в, отримаемо:

д2

и

дх2

к+1

3=0,N

к+1

3=0,N

+ / (**+1)

Враховуючи апрокеимащю (5) заиишемо (14) у впглядк

¿(£+1)

а

/ (**+1)

д2и к+1 1 ик

дх2 0 О2 т

д2и к+1 1 ик+1 - и дг и%

дх2 N О2 т

+

+

Шдетавпмо (15) 1 (16) у (11) 1 (12), гЛ/шогЛдно:

2 +0(Т),

+ 0(т).

ди дх

к+1 О

ик+1 - ик+1 ¡1к+1 ик+1 - 4 / (гк+1) ¡1к+1

кк+1

2 о2

т

2 о2

+ 0

У +

т

(2.7)

(2.8)

(2.9) (2.10)

(2.11) (2.12)

(2.13)

ди дх

k+l .,k+1 „,fc+1 i fe+i „,fc+l „,fc f í+k+l\ uk+l

JV

/ife+1 2o2 r 2o2

Таким чином, отримали апрокеимацш другого порядку i будемо викориетову-вати 'íi*, а не апрокеимацш (8) першого порядку, Зазначимо тшьки, що шар (j = 0) вщповщае . lir.iii (шлыпп) меж i облип i (xq = sfe), a (j = N) - пра г,i ¡i (.r y = L), 11 i. t-етавимо (17) y (2), отримаемо для граничних г,уз. lir, cítkh:

„MI _ 7,fc+1 hk+1 7,fc+1 _ „.k f (fk+l) hk+l

- - 2J + № = (ít+1).

Шеля перетворень:

fahM a A t+l ( -a \ M rf« , ,t+n f(tM)hM

(йГ + TF* - 3) + [ñ^) = - ** ) - « w ' a6°

+ ^r - a5+1 + ffl) aí+1 = - № (í'+1) ir -' ^ht+K

Введемо позначення:

2o2 hk+l 2o2/? ^2o2 O0 = O, o0 = 7T-T H----, Со

hk+l т a ' /ife+1 '

do = —4 - (fio (tk+l) — - f (tk+l) hk+\ (2.15) т a

тощ

h UQ+1 + cqu\+1 = 4, i = 0. (2.16) Тепер шдетавимо (18) y (3) i впконаемо ана. юпчш перетворення:

[hk+1J uN-i + [hk+i + 2o2t +ó)un )+ 2aV UN 7 2a2 ,

4+Л + (Ц + ^ + 4+l = f ^ r1) + ^4-/r1)

Введемо позначення:

^2o2 , 2o2 hk+l 2a¿S

°'N ~ Uk+1 ' ^ ~~ 7T7T ^^---1--> CN — o,

¡l ri т 7

div = —^L (ífe+1) + —- / (ífe+1) (2.17)

7 4 ' т 3 урахуванням позначень, отримаемо:

о^+Д + = djv, j = ЛГ. (2.18)

Таким чином, маемо апрокеимащю граничних умов третього роду на вшьнш 1 правш границях з тим же порядком, що й апрокеимащя (7) диференцшного р1вняння (1), Перепишемо (7) у виглядк

(7^) и1+1 +(~ + 7^2 ) +1 + (7^) «?й = - + / ■

' и 3 \(Ь+1)У 3+1 г ^ >

I Гн. 1я введения позначень, будем о мати:

_„-' 1 О02 _ 2 ук

а3 = —Ь3 = ± + , с, = —43 = -А + / (2.19)

'' (кк+1Т 1 г ' (Нк+1у ' (кк+1У ' г

ази^-1 + ьзи)+1 + сзи)Х1 = ¿3,3 = 1,^-1 (2.20)

Оетаточно, отримаемо шукану неявну р*1 шппеиу схему:

и® = г/; , ] = 0, Ж, (початкова умова);

Ъ0ик+1 + сом^1 = <1о,э = 0, (вшьна граииця); аIик ' | + 1>)ик ' 1 + = с^-,= 1, ЛГ — 1, (внутриншвузлп);

оагм^^Д + = = -/V", (права границя);

00Нт = {хк = 8к+з- Нк, Ьк = {!.- / = 5 (¿*) } ,

j = 0, Ж; Г = к-т, к = 0, М:

2 о2 2а2/? ^2о2

. 2а

2

"и0

о

т у ' а у '

_ ^2о2 2о2 2а28

Ч V — , ь , 1 , УЛГ — ТТТТ ^---'--' Слг —

Г1 Г1 т 7

Маемо систему л'нппннх алгебраТчних рпшяпь з трьох д1агональною матрицею, яка може бути розв'язана, наприклад, методом прогонки.

3. програмна реал13ац1я р13ницев01 схеми

Отриману кшдево-р1зницеву схему реализуемо засобами нрограмуван-пя MathCad (рис.1). Результатом обчиелень буде двом1рний масив Z, який вццювщае значениям шукано! функди ик = и (xk,tk) у вузлах еггки.

г -

for ie 0.. N

к

for к е 0..N - 1

t< т-{к +1) b^t^W]

(aot-o ед (-о)

1А! h 2-A2f3(t) -2-А2 b-Uo,k ф0р)1А2

bn »--;— +---С()<--7— dftg-----F(t)-ll

h T q b т a

—2-a~ 3-az h 2-а:-й 2-az-tbL(t) ь-им.ь

- b\j —■--\---h- ---+--f(t)-b

b h т к ' -4 т

a\

for jcl_N I

J '

-A" , 1 2-A' — bj <- - + —r h2 1 T bz

-A Ui к

for ieo..N for j eo..N Rj j f- 0 (r.D;(i<-bg Sj,|t-ti rn.k 1< for ре 1 ..N - !

(Rp.p-l^ap Rp,p<-bp Rp,pi-I cp ) F v- lsoIve{R,d) for pe 0..N

Up rk+i Fp return U

Рис. 1. Програма рпиицево!' неявно! схеми (MathCad)

Зауважимо, що реа,;пзад1я неявно! схеми, обумовила викориетання квадратно!' ciTKH (N = М), що не внливае на зб1жшеть та етшюеть алгоритму. Збыьшення кыь-Kocri крок1в еггки збыьшуе точшеть отримуваного результату. Надамо пояснения ви-користаним (рис.1) нозначенням: N- число кроюв еггки но нроеторовш координат! (довжиш); М- число кроюв еггки но чаеовш координат! (N = М); L межа обласп но нроеторовш координата (х = L); Т- максима.1!ьне значения часу (t <Т < ос): А- ко-ефщ1еит а в piBimuui дифузи (визначаеться влаетивоетями середовища); т = Т/М-ноетшний часовий крок еггки; / (¿)- функд1я право!' частили р1вняния (1); S (¿)-функд1я выыю! ,;пво! межц ф (х)- функд1я ночатково! умови (4); фО (¿)- функд1я Ра (t) нраво! частили умови на выьшй мож1 (2); фЬ (t) - функд1я ipi (t) нраво! части-пи умови на нравш мож1 (3); а, ¡3- вццювццп коефщ1енти граничних умов (2); 7,5-вццювццп коефщ1енти граничних умов (3).

4, Приклад застосування програми

У якосп приклада, розв'язуватимемо паетуппу задачу:

д2и ди п , , — - — = 0, а (¿) < х < 1 , I > 0;

^ = ± (щ - и) * = « 00 = 0.5 - I2/8, ^ > 0, щ = 1; ди (1 I)

= 0. X = 1.1 > 0; и (х. 0) = 0. а (0) <х <1.1 = 0.

ах

У вццювццюсп з нозиачеииями у программ дня цього прикладу матимемо:

/(¿) = 0,Ь = 1,Г = 2,.4 = 1,^ (.т) = 0, 5 (¿) = 0.5 - ¿2/8, а = 1, ¡3 = ¡3 (¿) = г¿5 ОО/гЙ, ф0 (¿) = (¡Б ({)/&, т = = о,фЬ (¿) = 0.

На рис.2 наведено результата розрахупюв з урахуваппям порма.;пзади зпачепь дня проеторово! сггки. С,;пд зауважити, що можливосп еиетеми MathCad дозволяють нобудувати комбшацпо лшш р1впя та проетрово! сггки, що шдвищуе паглядшеть иодаппя отримапих результатов, Кр1м паведепих залежпоетей, вбудоваш можливосп MathCad дозво.няють нобудувати вццювццп поверхш та фупкдюпальш залежпоетг Програма дозволяе швидко змшювати фупкдпо, що визпачае характер выыю! межь отримувати вццювццп розподЬш та коптролювати зб1жшеть отримапого розв'язку. Наведена па рис.1 програма, може бути реа.;пзовапа будь якою мовою програмуваппя, за умови програмпо! реа.;пзади фупкди розв'язаппя еиетеми р1впяпь.

Рис. 2. Лши р1впя па р1зпидевш еггд1

Висновки

Кшцево-р1зницеву схему для ¡ндн'п Стефана побудовано так, що порядок апрокеимадп рпшяння збер1гаетьея i для апрокеимадп граничних умов, Такий iiiixii не тшькп впр1внюе та шдвищуе порядок апрокеимадп (однорщна апроксимащя), але й робить кшдево-р1зпидеву схему консервативною. Таким чином, враховуютьея vci види енергп, що викориетовувалися шд чае отримання рпшяння. Побудований алгоритм, дозволяе знаходити розв'язок для pi ¡них функдш, що впзначають вшьну гранпдю. Результата обчиелень легко r,iзуа. iiзуг.а ч п для будь-яких залежиоетей, За-зпачимо додатково, що викладений метод може бути заетосований до рпшяння (1) за наявноета конвективних члешв (пропорцшних noxi. niiii ди/дх ) та ч. lenir, проиор-щйних шукан'пЧ функдп.

Список литературы

1. В.А. Кудинов, Э.М. Карташов, В.В. Калашников. Аналитические решения задач тепломассопе-реноса и термоупругости для многослойных конструкций: учебное пособие.- М.: Высшая школа, 2005.- 429с.

2. Гурский Д.А., Турбина Е.С. Вычисления в Mathcad 12. — СПб.: Питер, 2006. — 544 с.

3. Формалев В.Ф., Ревизников Д.Л. Численные методы. - М.: ФИЗМАТЛИТ, 2004. - 400с.

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