http___blog.simmakers.com_quasilinear-heat-equation__lang_bg_1.xml
http___blog.simmakers.com_quasilinear-heat-equation_.xml
1.1 2.1 2.2 2.3 3.1 4.1 5.1 5.2 6.1 7.1 8.1 9.1 10.1 10.2 10.3 11.1 11.2 11.3 11.4 12.1 12.2 12.3 13.1 14.1 15.1 16.1 16.2 16.3 16.4 17.1 17.2 18.1 19.1 19.2 19.3 20.1 20.2 21.1 22.1 23.1 23.2 23.3 23.4 23.5 23.6 23.7 24.1 24.2 25.1 25.2 25.3 25.4 26.1 27.1 27.2 28.1 29.1 29.2 30.1 31.1 31.2 32.1 33.1 34.1 34.2 35.1 36.1 36.2 37.1 37.2 37.3 38.1 39.1 39.2 39.3 39.4 39.5 39.6 39.7 40.1 41.1 41.2 41.3 41.4 41.5 42.1 43.1 44.1 45.1 45.2 45.3 46.1 46.2 46.3 47.1 48.1 49.1 49.2 50.1 51.1 52.1 53.1 54.1 55.1 55.2 56.1 56.2 56.3 57.1 58.1 58.2 58.3 59.1 60.1 61.1 62.1 62.2 62.3 63.1 64.1 65.1 66.1 66.2 67.1 68.1 69.1 70.1 71.1 72.1 72.2 73.1 74.1 75.1 76.1 76.2 77.1 78.1 79.1 80.1 81.1 82.1 83.1 83.2 83.3 84.1 85.1 86.1 87.1 88.1 89.1 89.2 90.1 91.1 91.2 91.3 91.4 91.5 91.6 92.1 93.1 94.1 95.1 95.2 95.3 95.4 95.5 95.6 95.7 96.1 96.2 96.3 96.4 96.5 97.1 97.2 97.3 98.1 98.2 99.1 99.2 100.1 101.1 101.2 101.3 102.1 103.1 104.1 105.1 105.2 106.1 106.2 106.3 107.1 107.2 108.1 108.2 108.3 109.1 109.2 110.1 110.2 110.3 110.4 111.1 111.2 112.1 112.2 112.3 112.4 112.5 113.1 113.2 113.3 113.4 113.5 114.1 115.1 115.2 116.1 116.2 117.1 118.1 118.2 118.3 119.1 119.2 120.1 120.2 120.3 120.4
Можете да прочетете на оригиналната статия на уеб сайта на World Congress Engineering. Абстрактен — The quasilinear heat equation with thermal conductivity and heat capacity depending on the temperature field in three spatial dimensions is studied in application to the phase transition problem in permafrost ground. Условията, при които на редуващи се посоки Дъглас - Rachford краен разликата схема запазва числения стабилност са изрично формулирани. Сравнението с известния аналитично решение сходство на проблема Stefan в една пространствено измерение се извършва. Индекс Условия-Quasilinear уравнение топлина, Стефан проблем, крайни разлики, редуващи се посоки схема, числения стабилност. ВЪВЕДЕНИЕ От началото на формулировки на редуващи се посоки имплицитни методи (АДИ методи) [1], [2], те са изключително развита и е установено, голям брой приложения [3], [4]. Въпреки това, сериозни трудности се срещат с използването на тези методи при прилагане на проблемите със сложни геометрични [5] и / или нелинейни уравнения на математическата физика [6]. Въпреки, че на схемите на методите ADI се доказа, че е ефикасен и икономичен по отношение на разход на време и, в повечето случаи, безусловно стабилна, те имат някои недостатъци: Техните крайни разлика формулировки позволяват да се разгледа само правоъгълни пространствени домейни (поради commutativity изисквания, наложени на factorized и разделят оператори) [7]; Прилагането на ДДД схеми на проблемите с Neumann и Робин гранични условия, които са различни по време срещне сериозни проблеми, дължащи се на необходимостта от оценка или сближаване на тези гранични условия в междинните етапи на схемата [8]; Когато се използва за решаване на нелинейни уравнения топлина, операторите, представляващи схема ADI не пътуват, което води до загуба на безусловно стабилност на схемата [6]. Първият от посочените недостатъци могат да бъдат преодолени чрез използването на крайни методи елементи във връзка с техники за разделяне на оператора, или чрез техники на домейни на разпадане. Вторият и третият недостатъци представи важен проблем за успешното прилагане на схемата ADI. Доколкото ни е известно, не пълен анализ стабилност за схема за ADI прилага за нелинейно уравнение топлина в триизмерен пространствената област е наличен в литературата, като по този начин мотивира тази работа. Друга мотивация за настоящата работа е прилагането на схемата ADI за моделиране на топлопреминаване в мащабни екологични системи (e.g., large areas of permafrost ground) който, в случай на чисто изрични крайни разлики схеми, налага твърди ограничения на времето стъпка, за да се гарантира стабилността на цифровата. В същото време, прилагането на имплицитни схеми може да доведе до много по-големи разходи, отколкото тази на компютърни изрични схеми, особено за проблеми с бързо променящите се коефициенти в сложна геометрия и значително нехомогенни сита. По този начин, в моделиране на топлопреминаване в широкомащабни системи за необходимостта от полагането на оптимален избор между открити и тайни схеми възниква. В случай на крайни елементи метод, applied to modeling of processes in permafrost ground, анализ на цифровата стабилност изглежда толкова сложна, че стабилността критерий е често се създават емпирично [9]. In the present paper we discuss the application of the ADI Douglas – Rachford scheme to the solution of Stefan problem in porous permafrost ground. Изследването е организирано както следва: Следващият раздел съдържа проблема състав и някои предположения, които ще бъдат използвани за доказване на цифровата стабилност на схемата ADI, докато точка 3 излага себе си доказателство. Раздел 4 представя някои числени резултати и е последвано от заключенията. ПРОБЛЕМ настройка и формализъм Quasilinear топлина уравнение Уравнението описание за пренос на топлина в системата с фазов преход има следния вид: където е неизвестната функция (температура: при всяка фиксирана че е истинско линейно пространство на положителни стойности функции), е топлинен капацитет, е топлопроводимост. От ще означаваме с ограничена, открита подмножество на евклидовата пространство с граница , приключването на , означена . След класификацията в [10], ние наричаме Eq.(1) а quasilinear уравнение топлина. Началните и гранични условия се приема, че е: Считаме, в случаите, когато и , по този начин Eq. (1) е равномерно параболична и има уникален разтвор на началната-границата стойност проблем екв. Stefan problem in permafrost ground Базирайки се на резултатите от [12], ние формулираме модел за проблем Stefan без изрично позоваване на предната проследяване състоянието. This approach is justified by the observations that in case of explicit front-tracking models applied to Stefan problem in salted permafrost ground, там се появява един overcooled региона (замразени ресни зона) [13] (анализа на тази зона и свързаните замръзване възвишаем и криогенни процеси смукателни [9] са извън обхвата на настоящата работа). По този начин, вземем предвид фазата на прехода чрез въвеждане на ефективни топлинен капацитет който включва латентна топлина на единица маса където и are the values of volumetric heat capacities of the ground in thawed and frozen phases respectively, е част от замръзнала вода, на изглаждане параметър, е температурата на фазовия преход, е плътността на водата. Коефициентът на топлопроводност се приема за където и са стойностите на топлопроводнос в размразени и замразени фази съответно. Линейното пространство За формулирането на крайните разликата схема ще се въведе следната процедура дискретизация. Нека вектор има положителни координати и е съвкупност от всички точки , индексите са цели числа. Две точки и принадлежащи към се наричат съседите, ако . Точките всички чиито съсед принадлежат са обозначени с . Следващ [14], ще означаваме точките с имоти, които най-малко един съсед принадлежи към от . По този начин, пълния пространствен окото е , е набор от гранични точки (извън ). За следващия приемем, че е хомогенна (т.е., са константи) кубичен домейн. Нека и . След това време, окото се определя, както следва: Приблизителните решения на екв.(1), , са определени и да вземат своите стойности в реално ограничен тримерно линейно пространство , размерът на които е равен на броя на точките, в . A функцията се нарича допустима, ако и на . Както бе отбелязано в [14], е, изобщо, не е известен и по този начин следната хипотеза следва да се: съществува нулева последователност (че клони към нула) на окото разстояния, така че и винаги принадлежи на последователността . Това предположение е гаранция, че допустима функция е еднозначно определено в по отношение на началните и граничните стойности на . За всеки дефинираме първи ред напред и назад разлика оператори: където е броят на точките в заедно -посока в декартови координати. За всяко подходящо дефинирана функция ние се По този начин, са стойностите на при фиктивни междинни възли на мрежата в . От уравнения. (8)-(9) , се получава: където и съответства на стойността на в и се определя като . Ние също така предоставя с -вътрешен продукт и индуцирана норма на : Максималната норма се определя, както следва: Ние дефинираме на модифициран -вътрешен продукт и индуцирана норма от където и съответстват на стойностите на и в съответно. Внимание, че така определения оператор предполага, че обикновено не е вярно, че (виж дискусия в следващия подраздел). Редуването посока схема краен разликата Нека сега да постави в кореспонденция до Eq. (1) следните крайни разлики схема: където се дават от уравнение. Като се има предвид Eq. (16) може да бъде пренаписана в следния вид: където (еднакви за ), и са допълнителни функции, които се използват в междинно съединение steps.Note че индексите в първото уравнение на системата (17) промени в следните диапазони , , докато за последните две уравнения (17) , , и , , съответно. Схемата (17) е описано в [1] и се нарича Douglas-Rachford схема. Това беше посочено от Дъглас в [1] че схемата (17) с определена от уравнение. (11) на местно равнище е от втори ред, правилно в пространството, освен ако гранични условия понижаване на точността. От Eq. (16) може да се види, че Eq. (1) е приблизително със точност в . Границата и начални условия (3)-(2) са формулирани на окото по следния начин: където . Добре известно е, че схемата (16) е безусловно стабилна, когато и в уравнение. (15) са константи (т.е., и ) и са крайни разлика форми на Лаплас оператори. Това лесно се вижда от уравнение.(16) и Eq. (14) ако забелязали, че са selfadjoint положителни определени оператори, . В този случай може да се използва следната теорема се оказа в [15]: Теорема 1: Ако операторите и в ур.(16) са самостоятелно прилежаща, не се променят във времето, и е положително определена, След това състоянието е необходимо и достатъчно за следващата оценка за провеждане: където е норма, предизвикана от -вътрешен продукт (14) в случай на постоянно и , е първоначалното състояние, и е цифровата разтвор, получен в -тата итерация. Наистина, ако и са константи, операторите определена от уравнение. (15) са по двойки комутативен и по този начин операторът връзка (19) се отнася и за всички малки . Връзка (20) означава, че схемата (16) числено е стабилна, когато и са константи. Въвеждането на и които зависят от прави оператори определена от (15) нелинейни, Т.е.. че по принцип не е вярно, че . Този факт го прави невъзможно да се докаже числен стабилност в една норма, както това се прави в екв.(20). За да продължим, ние леко се промени определението на оператор и да се направят следните допускания: Предположение 1: В схемата ADI (16) операторите действат от векторите по следния начин: Това предположение води до разглеждане на линейна схема краен разлика от тип (16) за термична уравнение с променливи коефициенти и , формите на които зависи от температурата областта изчислени при предишното ниво на време. Този вид на предположението, обикновено се използва, когато се обмисля нелинейни зависими от времето проблеми (e.g., [16]). Нека да пренапише (21) в следната еквивалентна форма: където Предположение 2: Функциите , и (определена от (4)-(6), (9) и (24)) са набелязани елементи в последователност пространство , елементите на които отговарят на следните отношения: когато операторите са определени по същия начин да (8); са реални положителни константи. е Banach пространство с максимална норма. III. NUMERICAL СТАБИЛНОСТ НА СХЕМА ADI прилага към нелинейни УРАВНЕНИЕ HEAT В този раздел ще се възползват от Предположения 1 и 2 и да получи условията, при които схемата (16) с Eq.(21) може да се прилага успешно за решаване на нелинейни топлина уравнение (1). Основният проблем, трябва да се разгледа, е следното: получи някаква връзка между енергийните норми на операторите които биха могли да се гарантира стабилността на цифровата (16). За да продължим, нека приеме следното определение на числения стабилност, дадена в [15]: Определение 1: Нека е положително определена самостоятелно прилежаща оператор, който действа в (истински) Хилберт пространство H. Нека да бъде постоянен, така че за всеки от (7), е константа, не зависи от . Схемата (16) се казва, да бъде стабилна в пространството Хилберт (с скаларно произведение и норма, определена от (14)) Ако по някаква и за всеки достатъчно малък и , за решаване на уравнението следната връзка притежава: където е идентичността оператор. За да докаже, че схемата (16) е стабилен в горния смисъл, ние ще се нуждаят от някои помощни резултати доказани в [15]: Лема 1: Нека бъде Хилберт пространство (реално или комплекс). Нека е линеен оператор , така че . Тогава връзката (с ) е достатъчно да държи. Лема 2: Нека да бъде в зависимост от времето (за от (7)), положително определена и самостоятелно прилежаща оператор. След това, връзка е еквивалентна на връзка където е константа, която не зависи от . Сега, въз основа на резултатите от [15], можем да формулираме основната теорема: Теорема 2: Нека операторите (определена от (21))отговаря на следното условие за всеки където е постоянна. Тогава връзката с от (7), ) е достатъчно за следващата оценка за провеждане: където са разтвори на екв.(16). Доказателство: Операторите (определена от (21)) са положителни определено и самостоятелно прилежаща, по този начин позволява да се определи и . Нека да пренапише Eq.(16) в следния вид: където Сега, използвайки Лема 1 и (34) ние получаваме: По този начин, за и човек трябва: Като се има предвид Лема 2 и първото уравнение на (38), ние получаваме: Следователно, И, накрая, връщане на скаларни продукти, , се получава (35). Нека сега да разгледа условията, при които връзка (34) притежава. Човек може да види, че l.h.s на (34) съдържа условията и (освен ) които са необходими за неравенството да държи. Като се има предвид Успение Богородично 2, обикновени изчисления на количество ни дава: където на макс() и мин.() поето всички . От (44)-(48) следва, че за неравенството (34) да държи, на достатъчно условие е: или, по-строго: От (50), на достатъчен критерий стабилност (, която се отнася и и гарантира, че (34) държи всеки път, когато нивото ) може директно да бъде изразена както следва: където е константа, която може да бъде извлечена от (50). Имайте предвид, че (51) прилича на учебник израз за стабилност критерий на крайно-разлики схема за onedimensional уравнение топлина с променливи коефициенти:. Числени експерименти В този раздел ще покаже резултатите от прилагането на схемата ADI (16) за решаването на проблем Stefan в правоъгълна паралелепипед. Ние сравняваме получените числени резултати с известните аналитични изрази на приликата решения за безкрайност домейни [17]. Да разгледаме следния проблем: едномерен semiinfinite домейн съдържа материал в течно състояние в и в неговото твърдо състояние в , е точката, която разделя течността и твърдите фази. Температурата в момента в момента момент е обозначен с . Температурата на фазов преход се приема за градуса по Целзий. Ние също така да поемат нулеви топлинни потоци в краищата на и piecewise постоянна първоначално състояние: Формулираното проблем е известно, че имат следния разтвор [17]: където и постоянната се получава чрез решаване на уравнението: където ERFC() е допълваща функция грешка, е латентна топлина на фазов преход, и са топлинните diffusivities в твърдата и течната фаза съответно които се приема за константа във времето и позицията. В цифровата нашия експеримент сме счита пространствената област с размери: , м. Пространственото окото има следните номера възли: 4х4х172 възли. На всички лица на паралелепипеда условия на гранични Neumann с нулеви топлинни потоци са наложени. Първоначалното състояние се приема, че е (52) с м. Топлопроводнос и топлинни мощности на размразени и замразени фази, които се появяват в (4)-(6) се приема, че е: , , . Числено стабилност за проблем с тези параметри е постигната, когато времето стъпка и пространствено стъпка са били обект на следното условие: . Това ни даде броя на повторения: 6486. Смокиня. 1. В зависимост от температурата на Z-координата: Сравнението на числовите резултати (сини звезди и зелени триъгълници), получени със схемата (16), и аналитичен израз на сходството разтвор (53) (червени кръгове). Време момент 28 дни. Вградени показва същите криви при температури под TPH = 0 Целзий. Red вертикална линия от вътрешната страна на паралелепипед показва набор от точки, , за което резултатите са представени. Смокиня. 1 показва числовите резултати по линията м (т.е., по протежение на централната линия на популярното през паралелепипед) заедно със съответната аналитичния резултат, получен с (53). Човек може да види както и качествен и количествен споразумение между числени и аналитични резултати, точността е в рамките на 0.5 градуса по Целзий. Както се очакваше, по-добро споразумение се постига при по-големи стойности на параметър изглаждане в (5), който е отговорен за сметка на остротата на фазовия преход. Един добър споразумение между числени и аналитични резултати се наблюдава и за времето, еволюцията на температурата в дадена точка. По този начин, Смокиня. 2 показва времето зависимост от температурата в точката с координати м. Смокиня. 2. Зависимостта на температурата от времето тон в точката м: Сравнението на числовите резултати, получени със схемата (16), (сини звезди и зелени триъгълници) и аналитичен израз на сходството разтвор (53) (червени кръгове). дни. ИЗВОДИ Ние показахме, че при допусканията 1 и 2, единния цифров стабилност (Определение 1) може да се установи за Douglas - Rachford схема (17) (която се доближава уравнението (1)), общата форма на достатъчен критерий стабилност е (51). Използваният подход се състои в редукция на quasilinear уравнението на уравнение с променливи коефициенти (при спазване на ограниченията (25)-(27)), формите на които се променят след всяка итерация. Сравнението на получените числени резултати и аналитично решение (Фиг.(1) - (2)) suggests that the ADI method can be successfully applied to the solution of Stefan problems in permafrost ground. БЛАГОДАРНОСТИ ИЗТОЧНИЦИ Douglas, "Алтернативна методи насока за три космически променливи,"Числени математика, полет. 4, не. 1, п.п.. 41-63, Декември. 1962. Yanenko, Методът на дробни стъпки;: Решаването на проблемите на математическата физика в няколко променливи, 1во изд. Springer-Verlag, 1971. Марчук, Разделяне и променлив Посока Методи, бъде. Краен Методи Разлика. Северна Холандия: Elsevier Science Publishers B.V., 1990, полет. 1, п.п.. 199-462. J. в "тон Hout и B. Welfert, "Безусловната стабилност на схемите secondorder ADI прилага за многомерни дифузионни уравнения със смесени производни условия,"Приложна Числени математика, полет. 59, не. 3-4, п.п.. 677-692, Море. 2009. Ангот, J. Keating, и Р. Минев, "A алгоритъм посока разделяне за несвиваем поток в сложна геометрия,"Компютърни методи в приложна механика и инженерство, полет. 217-220, п.п.. 111-120, 2012. Hansen и A. Остерман, "Dimension разделяне за quasilinear параболични уравнения,"IMA вестник Числени методи, полет. 30, не. 3, п.п.. 857-869, 2010. Бъркоф и R. Varga, "Косвени редуващи методи посока,"Trans. Amer. Math. Białecki и R. Fernandes, "Променливо-посока имплицитно правоъгълната схема сплайн колокиране за нелинейни параболични проблеми на правоъгълни полигони,"SIAM вестник за научно Computing, полет. 28, не. 3, п.п.. 1054-1077, 2006. Thomas, P. Cleall, И. Li, C. Harris, и М. Основна Lütschg, “Modelling of cryogenic processes in permafrost and seasonally frozen ground,"G'eotechnique, полет. 59, не. 3, п.п.. 173-184, 2009. Ladyzenskaja, В. Solonnikov, и N. Ural'ceva, Линейни и квази-линейни уравнения от параболичен тип, бъде. Преводи на Математическа Монографии. Провидение, RI: Американския Математическо общество, 1968, полет. 23. Kamenomostskaya, "По въпроса Стефан." Mat. Samarskii and P. Vabishchevich, Математическо моделиране, Обем 1, Изчислителна Heat Transfer, обем 1 и. Bronfenbrener и E. Korin, "Двуфазен условия зона образуване по замразяване на порести среди,"Вестник на кристалния растеж, полет. 198-199, п.п.. 89-95, 1999. Чета, "Алтернативна посока и полу-ясна разлика методи за параболични частни диференциални уравнения,"Числени математика, полет. 3, не. 1, п.п.. 398-412, 1961. Samarskii and A. Gulin, Стабилност на Разлика Схеми. Moskow: Наука, 1973. Денди, "Алтернативна методи посока нелинейни Time- Зависимите проблеми,"SIAM J. Номер. Javierre, C. Vuik, F. Vermolen, и S. Ван дер Zwaag, "Сравнението на числени модели за едномерни проблеми Стефан,"Вестник на компютърна и приложна математика, полет. 192, не. 2, п.п.. 445-459, 2006.
1.1 2.1 2.2 2.3 3.1 4.1 5.1 5.2 6.1 7.1 8.1 9.1 10.1 10.2 10.3 11.1 11.2 11.3 11.4 12.1 12.2 12.3 13.1 14.1 15.1 16.1 16.2 16.3 16.4 17.1 18.1 19.1 19.2 19.3 20.1 20.2 21.1 22.1 23.1 23.2 23.3 23.4 23.5 23.6 23.7 24.1 24.2 25.1 25.2 25.3 25.4 26.1 27.1 27.2 28.1 29.1 30.1 31.1 31.2 32.1 33.1 34.1 34.2 35.1 36.1 37.1 37.2 37.3 37.4 38.1 38.2 39.1 39.2 39.3 39.4 39.5 39.6 40.1 41.1 41.2 41.3 41.4 41.5 42.1 43.1 44.1 45.1 45.2 45.3 46.1 46.2 46.3 47.1 48.1 49.1 49.2 50.1 51.1 52.1 53.1 54.1 54.2 55.1 55.2 56.1 56.2 56.3 57.1 58.1 58.2 59.1 60.1 61.1 62.1 62.2 62.3 63.1 64.1 65.1 66.1 66.2 67.1 68.1 69.1 70.1 71.1 72.1 72.2 73.1 74.1 75.1 76.1 76.2 77.1 78.1 79.1 80.1 81.1 82.1 83.1 83.2 83.3 84.1 85.1 86.1 87.1 88.1 89.1 89.2 90.1 91.1 91.2 91.3 91.4 91.5 91.6 92.1 93.1 94.1 95.1 95.2 95.3 95.4 95.5 95.6 95.7 96.1 96.2 96.3 96.4 97.1 97.2 97.3 98.1 98.2 99.1 99.2 100.1 101.1 101.2 101.3 102.1 103.1 104.1 105.1 105.2 106.1 106.2 106.3 107.1 108.1 109.1 110.1 110.2 110.3 110.4 111.1 112.1 113.1 113.2 113.3 114.1 114.2 114.3 115.1 116.1 117.1 118.1 118.2 119.1 119.2 120.1
You can read the original article at the World Congress Engineering web site. Abstract — The quasilinear heat equation with thermal conductivity and heat capacity depending on the temperature field in three spatial dimensions is studied in application to the phase transition problem in permafrost ground. The conditions under which the alternating directions Douglas – Rachford finite difference scheme retains numerical stability are explicitly formulated. The comparison with the known analytical similarity solution to the Stefan problem in one spatial dimension is performed. Index Terms—quasilinear heat equation, Stefan problem, finite differences, alternating directions scheme, numerical stability. INTRODUCTION Since the early formulations of the alternating directions implicit methods (ADI methods) [1], [2], they have been tremendously developed and found a vast number of applications [3], [4]. Nevertheless, serious difficulties are encountered with the use of these methods in application to problems with complex geometries [5] and/or nonlinear equations of mathematical physics [6]. Although the schemes of the ADI methods are proved to be efficient and economic with respect to time consumption and, in most cases, unconditionally stable, they exhibit some disadvantages: 1) Their finite-difference formulations permit to consider only rectangular spatial domains (due to commutativity requirements imposed on the factorized and split operators) [7]; 2) The application of the ADI schemes to the problems with Neumann and Robin boundary conditions that are varying in time encounters serious problems due to the necessity of evaluation or approximation of these boundary conditions at the intermediate steps of the scheme [8]; 3) When applied to the solution of nonlinear heat equations, the operators constituting an ADI scheme do not commute, thus leading to the loss of unconditional stability of the scheme [6]. The first of the above disadvantages can be overcome either by the use of finite elements methods in conjunction with operator splitting techniques, or by domain decomposition techniques. The second and the third disadvantages present an important problem for the successful application of the ADI scheme. To the best of our knowledge, no complete stability analysis for an ADI scheme applied to the nonlinear heat equation in a three-dimensional spatial domain is available in the literature, thus motivating this work. Another motivation for the present work is the application of ADI scheme to the modeling of heat transfer in large scale environmental systems (e.g., large areas of permafrost ground) which, in the case of purely explicit finite-differences schemes, imposes stiff constraints on the time-step in order to guarantee the numerical stability. At the same time, implementation of implicit schemes can often lead to much greater computing expenses than that of explicit schemes, especially for the problems with rapidly changing coefficients in complex geometries and substantially nonhomogeneous meshes. Thus, in modeling of heat transfer in large scale systems the necessity of making an optimal choice between explicit and implicit schemes arises. In case of finite-elements method, applied to modeling of processes in permafrost ground, the analysis of numerical stability appears to be so complex, that the stability criterion is often established empirically [9]. In the present paper we discuss the application of the ADI Douglas – Rachford scheme to the solution of Stefan problem in porous permafrost ground. The paper is organized as follows: next section contains the problem formulation and some assumptions that will be used in the proof of numerical stability of the ADI scheme while section 3 exposes the proof itself. Section 4 presents some numerical results and is followed by Conclusions. PROBLEM SETUP AND FORMALISM Quasilinear heat equation The equation describing heat transfer in a system with phase transition has the following form: where is the unknown function (temperature: at every fixed being a real linear space of positive valued functions), is the heat capacity, is thermal conductivity. By we denote a bounded, open subset of an Euclidian space with boundary , the closure of being denoted by . Following the classification given in [10], we call Eq.(1) a quasilinear heat equation. The initial and boundary conditions are taken to be: We consider the cases when and , thus Eq. (1) is uniformly parabolic and has a unique solution to the initial-boundary-value problem Eq. (3)–(2).[11]. Stefan problem in permafrost ground Basing on results of [12], we formulate the model for the Stefan problem without explicitly invoking the front-tracking condition. This approach is justified by the observations that in case of explicit front-tracking models applied to Stefan problem in salted permafrost ground, there appears an overcooled region (frozen fringe zone) [13] (the analysis of this zone and related frost heave and cryogenic suction processes [9] are beyond the scope of the present work). Thus, we take into account the phase transition by introducing effective heat capacity which incorporates the latent heat per unit mass where and are the values of volumetric heat capacities of the ground in thawed and frozen phases respectively, is the fraction of frozen water, is the smoothing parameter, is the phase transition temperature, is the density of water. The thermal conductivity is taken to be where and are the values of thermal conductivities in thawed and frozen phases respectively. The linear space For the formulation of finite difference scheme we introduce the following discretization procedure. Let the vector have positive coordinates and be the set of all points , the indices being integers. Two points and belonging to are called neighbors if . The points all of whose neighbor belong to are denoted by . Following [14], we denote the points with the property that at least one neighbor belongs to by . Thus, the full spatial mesh is , being the set of boundary points (outside of ). For the following we assume that is a homogeneous (i.e., are constants) cubic domain. Let and . Then the time mesh is defined as follows: The approximate solutions of Eq.(1), , are defined on and take their values in a real finite-dimensional linear space , the dimension of which is equal to the number of points in . A function is called admissible if and on . As was pointed out in [14], is, in general, not known in and thus the following assumption should be made: there exists a null sequence (that converges to zero) of mesh spacings such that and always belongs to the sequence . This assumption guarantees that an admissible function is uniquely specified in in terms of the initial and boundary values of . For every we define the first-order forward and backward difference operators: where being the number of points in along -direction in the Cartesian coordinates. For any suitably defined function we set Thus, are the values of at the fictitious intermediate nodes of the mesh at . From the Eqs. (8)–(9) one obtains: where and corresponds the value of at and being defined as . We also provide with the -inner product and induced norm on : The maximum norm is defined as following: We define on a modified -inner product and induced norm by where and correspond the values of and at respectively. Note, that the so defined operator implies that it is generally not true that (see a discussion in the next subsection). Alternating direction finite difference scheme Let us now put in correspondence to Eq. (1) the following finite differences scheme: where are given by Eq. (15), , and at . Taking into account Eq. (8)–(11) and (15), Eq. (16) can be rewritten in the following form: where (same for ), and are the auxiliary functions that are used at the intermediate steps. Note that the indices in the first equation of system (17) change in the following ranges , , while for the last two equations in (17) , , and , , respectively. The scheme (17) has been described in [1] and is called the Douglas-Rachford scheme. It has been pointed out by Douglas in [1] that the scheme (17) with defined by Eq. (11) is locally second-order correct in space, unless the boundary conditions lower the accuracy. From Eq. (16) one can see that Eq. (1) is approximated with the accuracy in . The boundary and initial conditions (3)–(2) are formulated on the mesh in the following way: where . It is well known that the scheme (16) is unconditionally stable when and in Eq. (15) are constants (i.e., and ) and are the finite-difference forms of the Laplace operators. This is readily seen from Eq.(16) and Eq. (14) if noticed that are selfadjoint positive definite operators, . In this case one can use the following Theorem proved in [15]: Theorem 1: If operators and in Eq.(16) are self-adjoint, do not change in time, and is positive definite, then the condition is necessary and sufficient for the following estimate to hold: where is the norm induced by the -inner product (14) in the case of constant and , is the initial condition, and is the numerical solution obtained at the -th iteration. Indeed, if and are constants, the operators defined by Eq. (15) are pairwise commutative and thus operator relation (19) holds for any small . Relation (20) means that the scheme (16) is numerically stable when and are constants. The introduction of and that depend on makes the operators defined by (15) nonlinear, i.e. it is generally not true that . This fact makes it impossible to prove numerical stability in a single norm, as it is done in Eq.(20). In order to proceed, we slightly modify the definition of operator and make the following assumptions: Assumption 1: In the ADI scheme (16) the operators are acting on the vectors in the following way: This assumption leads to consideration of a linear finite difference scheme of type (16) for the heat equation with variable coefficients and , the forms of which depend on the temperature field calculated at the previous time level. This kind of assumption is commonly used when considering nonlinear time dependent problems (e.g., [16]). Let us rewrite (21) in the following equivalent form: where Assumption 2: The functions , and (defined by (4)–(6), (9) and (24)) are mapping the elements into a sequence space , the elements of which satisfy the following relations: where the operators are defined similarly to (8); are real positive constants. is the Banach space with the maximum norm. III. NUMERICAL STABILITY OF ADI SCHEME APPLIED TO THE NONLINEAR HEAT EQUATION In this section we make use of Assumptions 1 and 2 and obtain conditions under which the scheme (16) with Eq.(21) can be successfully applied to the solution of nonlinear heat equation (1). The main problem one has to consider is the following: obtain some relation between energy norms of operators that could guarantee numerical stability of (16). In order to proceed, let us accept the following definition of numerical stability given in [15]: Definition 1: Let be a positive definite self-adjoint operator acting in a (real) Hilbert space H. Let be a constant, such that for any from (7), being a constant that does not depend on . The scheme (16) is said to be stable in a Hilbert space (with scalar product and norm defined by (14)) if for any and for any sufficiently small and , for the solution of equation the following relation holds: where is the identity operator. In order to prove that the scheme (16) is stable in the above sense, we will need some auxiliary results proved in [15]: Lemma 1: Let be a Hilbert space (real or complex). Let be a linear operator on , such that . Then the relation (with ) is sufficient for to hold. Lemma 2: Let be a time-dependent (for from (7)), positive definite and self-adjoint operator. Then, relation is equivalent to relation where is a constant that does not depend on . Now, based on the results of [15], we can formulate the main theorem: Theorem 2: Let the operators (defined by (21))satisfy the following condition for any where is a constant. Then the relation with from (7), ) is sufficient for the following estimate to hold: where are the solutions of Eq.(16). Proof: Operators (defined by (21)) are positive definite and self-adjoint, thus permitting to define and . Let us rewrite Eq.(16) in the following form: where Now, using Lemma 1 and (34) we obtain: Thus, for and one has: Taking into account Lemma 2 and the first equation of (38), we obtain: Hence, And, finally, returning to the scalar products, one obtains (35). Let us now examine the conditions under which relation (34) holds. One can see that the l.h.s of (34) contains the terms and (in addition to ) which are required for the inequality to hold. Taking into account Assumption 2, a straightforward calculation of a quantity gives us: where the max() and min() being taken over all . From (44)–(48) it follows that for the inequality (34) to hold, the sufficient condition is: or, more rigorously: From (50), the sufficient stability criterion (that relates and and guarantees that (34) holds on every time level ) can be directly expressed as follows: where is a constant that can be deduced from (50). Note that (51) resembles the textbook expression for the stability criterion of a finite-differences scheme for a onedimensional heat equation with variable coefficients:. NUMERICAL EXPERIMENTS In this section we demonstrate the results of application of the ADI scheme (16) to the solution of a Stefan problem in a rectangular parallelepiped. We compare the obtained numerical results with the known analytical expressions of the similarity solutions for infinite domains [17]. Consider the following problem: a one-dimensional semiinfinite domain contains the material in its liquid state in and in its solid state in , being the point that separates the liquid and the solid phases. The temperature at the point at time moment is denoted by . The phase transition temperature is assumed to be degrees Celsius. We also assume zero heat fluxes at the ends of and a piecewise constant initial condition: The formulated problem is known to have the following solution [17]: where and the constant is obtained by solving the equation: where erfc() is the complementary error function, is the latent heat of phase transition, and are the thermal diffusivities in the solid and the liquid phases respectively which are assumed to be constant in time and position. In our numerical experiment we have considered the spatial domain with dimensions: , meters. The spatial mesh has the following nodes numbers: 4x4x172 nodes. On all the faces of the parallelepiped the Neumann boundary conditions with zero heat fluxes have been imposed. The initial condition was taken to be (52) with meters. Thermal conductivities and heat capacities of thawed and frozen phases that appear in (4)–(6) are taken to be: , , . Numerical stability for the problem with these parameters has been achieved when the time step and spacial step were subject to the following condition: . This gave us the number of iterations: 6486. The dependence of temperature on the Z-coordinate: comparison of numerical results (blue stars and green triangles), obtained with the scheme (16), and analytical expression of the similarity solution (53) (red circles). Time moment 28 days. Inset shows the same curves at temperatures below Tph = 0 Celsius. Red vertical line inside the parallelepiped indicates the set of points, for which the results are presented. Fig. 1 shows the numerical results along the line meters (i.e., along the central line ofthe parallelepiped) together with the corresponding analytical result obtained with (53). One can see both a qualitative and quantitative agreement between numerical and analytical results, the accuracy being within 0.5 degree Celsius. As expected, a better agreement is achieved at greater values of the smoothing parameter in (5), which is responsible for the account of sharpness of the phase transition. A good agreement between numerical and analytical results is also observed for the time evolution of temperature at a given point. Thus, Fig. 2 shows the time dependence of temperature at the point with coordinates meters. The dependence of temperature on time t at the point meters: comparison of numerical results, obtained with the scheme (16), (blue stars and green triangles) and analytical expression of the similarity solution (53) (red circles). days. CONCLUSIONS We have shown that under the Assumptions 1 and 2, the uniform numerical stability (Definition 1) can be established for the Douglas – Rachford scheme (17) (which approximates the equation (1)), the general form of the sufficient stability criterion being (51). The used approach consists in reduction of a quasilinear equation to an equation with variable coefficients (subject to the constraints (25)–(27)), the forms of which change after each iteration. The comparison of obtained numerical results and an analytical solution (Fig.(1) - (2)) suggests that the ADI method can be successfully applied to the solution of Stefan problems in permafrost ground. ACKNOWLEDGMENT REFERENCES [1] J. Douglas, “Alternating direction methods for three space variables,” Numerische Mathematik, vol. [2] N. N. Yanenko, The method of fractional steps;: The solution of problems of mathematical physics in several variables, 1st ed. Springer-Verlag, 1971. [3] G. I. Marchuk, Splitting and Alternating Direction Methods, ser. Finite Difference Methods. North Holland: Elsevier Science Publishers B.V., 1990, vol. [4] K. J. in ’t Hout and B. D. Welfert, “Unconditional stability of secondorder ADI schemes applied to multi-dimensional diffusion equations with mixed derivative terms,” Applied Numerical Mathematics, vol. [5] P. Angot, J. Keating, and P. D. Minev, “A direction splitting algorithm for incompressible flow in complex geometries,” Computer Methods in Applied Mechanics and Engineering, vol. [6] E. Hansen and A. Ostermann, “Dimension splitting for quasilinear parabolic equations,” IMA Journal of Numerical Analysis, vol. [7] G. Birkhoff and R. S. Varga, “Implicit alternating direction methods,” Trans. Amer. Math. Soc., vol. [8] B. Bialecki and R. Fernandes, “An alternating-direction implicit orthogonal spline collocation scheme for nonlinear parabolic problems on rectangular polygons,” SIAM Journal on Scientific Computing, vol. [9] H. R. Thomas, P. Cleall, Y. C. Li, C. Harris, and M. Kern-Luetschg, “Modelling of cryogenic processes in permafrost and seasonally frozen ground,” G´eotechnique, vol. [10] O. A. Ladyˇzenskaja, V. A. Solonnikov, and N. N. Ural’ceva, Linear and quasi-linear equations of parabolic type, ser. Translations of Mathematical Monographs. Providence, RI: American Mathematical Society, 1968, vol. [11] S. Kamenomostskaya, “On the Stefan problem.” Mat. Sb., N. Ser., vol. [12] A. A. Samarskii and P. N. Vabishchevich, Mathematical Modelling, Volume 1, Computational Heat Transfer, volume 1 ed. [13] L. Bronfenbrener and E. Korin, “Two-phase zone formation conditions under freezing of porous media,” Journal of Crystal Growth, vol. [14] M. Lees, “Alternating direction and semi-explicit difference methods for parabolic partial differential equations,” Numerische Mathematik, vol. [15] A. A. Samarskii and A. V. Gulin, Stability of Difference Schemes. Moskow: Nauka, 1973. [16] J. E. Dendy, “Alternating direction methods for nonlinear Time- Dependent problems,” SIAM J. Numer. Anal., vol. [17] E. Javierre, C. Vuik, F. J. Vermolen, and S. van der Zwaag, “A comparison of numerical models for one-dimensional stefan problems,” Journal of Computational and Applied Mathematics, vol.