Численное моделирование теплообмена с фазовыми превращениями в задачах нефтегазопромысловой криологии тема автореферата и диссертации по физике, 01.04.14 ВАК РФ



На правах рукописи

СИГУНОВ Юрий Александрович

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

01.04.14 - теплофизика и молекулярная физика

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Тюмень - 1995

Работа выполнена на кафедре прикладной математики и автоматизированных систем управления Сургутского государственного

университета

Научный руководитель: член-корр. Инженерной Академии,

доктор технических наук, профессор Р.И.Медведский

Официальные оппоненты: член-корр. Академии Естественных наук,

доктор физико-математических наук, К.М.Федоров

кандидат технических наук, старший научный сотрудник Ю.С.Даниэлян

Ведущая организация: Тюменский государственный нефтегазо-

вый университет

Защита состоится " " СС НУП 'Э^Ъъ ¡X 1995 г. в 14 час. 30 мин. на заседании диссертационного совета К 064.23.02 в Тюменском государственном университете по адресу: 625003, г. Тюмень, ул. Семакова, 10. ИЧ

С диссертацией можно ознакомиться в библиотеке Тюменского государственного университета

Автореферат разослан " М " ^ H7JISL 1995 ,

Ученый секретарь диссертационного совета канд. физ.-мат. наук, с.н.с. Н.И.Куриленко

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы. Освоение нефтегазовых месторождений Западной Сибири поставило целый ряд проблем, от разрешения которых зависит эффективность строительства и надежность эксплуатации объектов добычи и транспортировки нефти и газа в условиях распространения вечномерзлых пород. Актуальность этих проблем возрастает в последнее время в связи с началом разработки новых крупных газоконденсатных месторождений п-ва1 Ямал, характеризующихся осложненными геокриологическими условиями. Важной составной частью возникающих задач является прогноз динамики и оценка последствий фазовых превращений воды и льда, происходящих в скважине и окружающих породах при их тепловом взаимодействии. Выбор технологии строительства и эксплуатации скважин на вечномерзлых грунтах требует предварительного многовариантного анализа теплообменных процессов, сопровождающихся фазовыми переходами, часто связанного с долгосрочным прогнозом, учитывающего большое число факторов, разнообразие строения и теплофизических параметров пород мерзлого разреза, различие конструкций и режима работы скважин. Удаленность объекта исследования на десятки и сотни метров от дневной поверхности существенно затрудняет прямые натурные наблюдения динамики фазовых превращений в системе скважина - породы. Все это повышает роль методов математического моделирования при решении обсуждаемого круга проблем.

Математическое описание теплообмена, сопровождающегося фазовыми превращениями, формулируется краевыми задачами с подвижными границами (задачами типа Стефана), представляющими собой сложный класс задач математической физики, имеющий разнообразные приложения в научных и технических дисциплинах. Возможности аналитического решения данного класса задач крайне ограничены, а применение численных методов наталкивается на серьезные трудности. Выделение положений фронтов фазового перехода, требующееся для корректного решения задачи, приводит к необходимости численного интегрирования уравнений в нерегулярных областях, разделенных движущимися границами локализации скрытой теплоты фазовых превращений, которые являются поверхностями сильного разрыва. Несмотря на большое число исследований методы численного решения подобного класса задач в настоящее время недостаточно разработаны.

Целью данной работы являлась разработка эффективных методов приближенно-аналитического и численного решения задач тепломассопереноса с фазовыми превращениями и применение разработанных методов к моделированию динамики фазовых переходов на примере задач, поставленных

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

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

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

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

Аппробаггия работы и публикации. Основные результаты докладывались и обсуждались на секции расширенного заседания Научного совета по криологии Земли АН СССР (Москва, 1991 г.); на Всесоюзном отраслевом совещании Газпрома СССР по вопросам разработки Бованенковского ГКМ (Москва, 1992 г.); на научных семинарах кафедры прикладной математики Сургутского государственного универитета под руководством профессора Г.И.Назина (Сургут, 1994-1995 г.г.); на научном семинаре по механике многофазных сред Тюменского научного центра СО РАН под руководством профессора А.А.Губайдуллина (Тюмень, 1995 г.).

По результатам диссертации опубликовано 10 печатных работ.

Структура и объем диссертапии. Диссертация состоит из введения, 6-ти глав, заключения и содержит 190 страниц сквозной нумерации, включая 9 таблиц, 39 рисунков, 116 библиографических ссыпок.

Автор выражает искреннюю признательность д.т.н., проф. Р.И.Медвед-скому, оказавшему определяющее влияние на направление и методологию выполненных в работе исследований.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

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

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

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

но на примере внутренней двухфазной задачи Стефана, представляющем наиболее сложный случай для аналитических подходов. Рассматривается затвердевание тела в форме плиты, цилиндра или сферы конечных размеров (О < £ < 1) с постоянной начальной температурой и граничным условием

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

Приближенное решение преследует цель определения координаты фазового фронта т]{т) и температуры поверхности тела ч/(х)=6\(\ьт) в произвольный момент времени. Уравнения для их определения выводятся на основе интегральных тождеств, выражающих баланс тепла в областях, занимаемых каждой из фаз вещества, и учитывающих энергию фазовых превращений. В работе получены следующие тождества

ат ах Кт а г

1

йт с1т Кт йг <>

V

в которых <р, - интегралы по пространственной координате, содержащие неизвестные распределения температур 9Х 0=0,1) и <92 (£,т) 0=2),

А(£) = (к = 0.1.2 соответственно для плоской, цилиндрической и сферической симметрии), \{г]) - весовая функция, используемая при интегрировании, вт - безразмерная температура фазового перехода. Безразмерные параметры определены выражениями

Р=ах/аг, = КМ=Л 1АТ/(а1 рЬ) ,

где а, и Х^ - температуро- и теплопроводность /-ой фазы, р Ь- теплота фазового перехода.

Переход от (1), (2) к уравнениям приближенного решения осуществляется с помощью подходящих аппроксимаций для температур 9У и в2 в интегралах <р0,<р1,<р2- Важно отметить, что тождества (1), (2) не содержат градиентов температуры, входящих в условия на фазовой границе и на возмущающей поверхности и исключаемых при выводе тождеств. Это

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

Дня построения приближений неизвестных температур использован предложенный U.C. Лейбензоном подход, когда аппроксимирующая функция находится из решения последовательности "близких" краевых задач в областях, занятых каждой фазой, при различных фиксированных положениях фазовой границы. Для вх (£,г) это соответствует принятою квазистационарного приближения

9l=0l = en+(V-9n)(\-V^)/V(rj)).

Распределение температуры 6г (£,г) в исходной фазе аппроксимируется решением соответствующей краевой задачи в области О <,Ç йт] с неподвижной границей щ при каждом фиксированном ее значении, играющем роль параметра. В связи с тем что точное решение этой задачи в случае конечной геометрии имеет сложное представление, в работе даны удобные в расчетах приближения для интегрального среднего от температуры по области 0 <, t] в виде функции времени, через которое выражается слагаемое pdç2l d г в (1), (2).

В конечном итоге приближенное решение задачи сводится к интегрированию системы двух обыкновенных дифференциальных уравнений относительно функций г](т) и ^т), описывающих движение фазовой границы и изменение температуры на возмущающей поверхности. Аналогичная система уравнений получается для внешней задачи Стефана. В общем случае их интегрирование производится численно. Для начального интервала времени получено асимптотическое приближение к решению.

Изложенная процедура вывода тождеств при определенных требованиях на аппроксимации температур приводит взамен (1), (2) к аналогичным ин-тегро-дифференциальным неравенствам. Основываясь на таком подходе в работе доказана теорема сравнения, обоснован способ получения двухсторонних оценок неизвестного закона движения фазового фронта и дана оценка погрешности предложенного метода, определяемая величиной безразмерных комплексов Km{\~0m) и Ктр вmsm. Разработанный метод аппробиро-ван на тестовых примерах в сравнении с решениями, полученными численными методами, и применен к решению ряда задач взаимодействия массива мерзлых пород со скважинами на различных этапах их функционирования.

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

При расчете динамики обратного промерзания пород вокруг остановленной скважины теплоотдача от нее моделировалась граничным условием Ш-го рода на ее контуре

£ = 1: ^=0^4 = 1,г). (3)

в котором безразмерный коэффициент теплоотдачи й(г) задавался различным в период работы скважины и при ее остановке. В период активного теплового воздействия работающей скважины температура жидкости в ней if/s(r ) задана. После остановки функция vs(r ) неизвестна и подлежит определению с учетом дополнительного уравнения

-c0di/fs/dT = h(rXWs-V), со = 0.5(суо)5/(ср),, (4)

являющегося следствием баланса тепла между остывающей скважиной и массивом пород.

На основе решения данной задачи проведено систематическое исследование скорости обратного промерзания вокруг остановленной скважины в зависимости от длительности предшествующего теплового воздействия, начальной температуры мерзлых пород, объемного содержания в них льда. Отмечается интенсивное поглощение тепла скважины и накопленного породами в окружающий массив с прекращением теплового воздействия, приводящее к началу обратного промерзания через 1-1,5 суток после окончания бурения. К этому моменту скважина теряет около 70% начального теплоза-паса. При прочих равных условиях величина времени возврата, отсчитываемая от момента остановки и отнесенная к длительности теплового воздействия составляет для различной начальной температуры пород: при -8°С ~ 1; при -5°С ~ 2; при -2°С - 7. Важным фактором, влияющим на сроки возврата фронта к контуру скважины, выступает льдистость пород. Общая закономерность, иллюстрируемая рис. I, заключается в уменьшении сроков возвра-

та фронта промерзания с понижением льдистости пород, несмотря на большую протяженность протаявшей зоны. Эта закономерность, впервые установленная нами по результатам вычислительных экспериментов, объясняет причину неравномерного замерзания ствола скважин, приводящего к замыканию водных объемов и опасному росту давления на колонны. Ранее А.В.Марамзиным указанная неравномерность связывалась лишь с наличием каверн, увеличивающих время промерзания оконоскважшшого пространства.

50

2 5

г.о

15

10

/ У"'

^ А -з

2

! \

5о ±,ат

Рис. I. Движение изотермы фазового перехода в породах с различным объемным содержание порового льда: 1 - 0.01; 2 - 0.05; 3 - 0.2; 4 - 0.4

Рис. 2. Движение границы фазового перехода в льдистом (/) и малольдистом (2) пластах при периодическом тепловом воздействии на породы

Важность предупреждения образования замерзающих замкнутых объемов обусловило рассмотреть в данной главе задачу о периодическом тепловом воздействии на породы, каковым может быть временная циркуляция в стволе скважины нагретой жидкостью. Прежде всего возникает вопрос о необходимой периодичности такого воздействия. С учетом отмеченного выше эффекта ее величина должна устанавливаться исходя из минимальной льдистости пород вдоль мерзлого разреза и планироваться на основании соответствующих расчетов. Расчетная модель остается прежней, а в условия (3), (4) закладываются последующие циклы теплового воздействия на породы с момента возврата изотерм фазовых превращений к контуру скважины в породах с наименьшей льдистостыо. На рис. 2 представлены результаты решения задачи для начальной температуры пород -2°С. Первый цикл воздействия соответствует бурению, последующие - прогреву скважины в течение 6 часов жидкостью нагретой до 60°С.

Картина криогенных процессов, вызванных остановкой скважины, дополняется далее исследованием динамики замерзания водосодержащих сред в окружении неограниченного массива мерзлых пород. Рассмотрены два таких примера: замерзание жидкости в стволе необсаженной скважины или

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

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

Условие на контакте скважины с породами выводится из баланса тепла: между породами и составным тепловым источником

~(ср)п AadTjdt- (ср\ AadTaJdt+AadQ/dt = ~2л Хх{гдТ lJdr)r,

между составными частями источника

-{ср\\d\ldt=2n

где Аа и Ац - площади, занятые в поперечном разрезе ствола скважины продавочной жидкостью и цементным кольцом, Та, \ и (ср)п, (ср)п - их температура и объемная теплоемкость, dQjdt - интенсивность тепловыделения при гидратации цемента.

Количество тепла, выделяющегося при гидратации цементного раствора, задавалось функциональной зависимостью (Запорожец ИД., Окороков С.Д., Парийский A.A.)

0(0 =1_

г

l + A2nJ/(rn)rff

л20

О

1

и-1

/(Г) = ехр^-^1п2, а + ЬТ

Smax

где О - максимальное выделяемое на единицу объема количество тепла.

cmix

остальные параметры носят экспериментальный характер

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

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

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

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

в котором функция А00 задает тип симметрии, а коэффициенты представляются кусочно-непрерывными функциями координаты и температуры. Кроме того, в любой точке х = z (г) справедлив закон сохранения энергии, выражение которого для рассматриваемого круга задач имеет вид

где величина pL(z,T) равна теплоте фазового перехода, если точка z(0 лежит на фазовом фронте, или полагается равной нулю в противном случае.

Для окрестности xa_t(t) < xn(t) < xn+1(r) произвольной точки области х = xn(t) вывод тождества осуществляется интегрированием уравнения (5) по областям jc„_, < х < хп и хи_, < х < хп соответственно с весовыми функциями А(х)ип_1(х) и А(х) vn(x), определенными равенствами

(5)

-(XdTjd х)хшг_0 + (Я д 7 ¡д х)х^ = pL(z,T)dz/dt,

(6)

х ли+)

/ч 1 с й а \ с йо . С (1 а

А» ■> Л(ст) Д. ' А(сг) Л А(сг)

Х„ ж лг„

и последующим суммированием результатов с учетом условия (6). Итоговое тождество имеет следующую форму

-ч т д 41 I

I (ср)—+ } (с/?) — V» {х)А{х)±с +Ап (р Ь)„ —^ = (8) * д\ 1 дх (1;

= Т" ГИг^"^ [ ¿ЗГ"*- А« = АМ (рЬ)я={рЫхя,ПхпУ).

Ап •> дх Д„_! * дх

Для граничных точек области аналогичные тождества выводятся интегрированием (5) по односторенней окрестности точки с учетом вместо (6) соответствующего граничного условия.

Использованием тождеств (8), форма которых не зависит от характера точки области, обходится необходимость включения в численную схему специальных уравнений, аппроксимирующих условия на фазовых границах. Это приводит к однородности и существенному упрощению вычислительного алгоритма, независимости его от числа подвижных фронтов и конкретной структуры расчетной сетки. Вследствие исключения прямых конечно-разностных аппроксимаций пространственных производных температуры на подвижных и возмущающих границах уменьшается погрешность решения, что позволяет существенно сократить необходимое число узлов сетки.

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

Температура на интервалах, ограниченных двумя последовательными узлами аппроксимируется квазистационарным приближением Т(х, 0 = уи {х) + у/п+1 ип (х) , хп < х < хпМ ,

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

а„ -*„-! + + ¿П + ) К + сп + ап + (Ьп + 2п ) + Сп =/п.х~/п,

в котором координата хп или температура у/п является известной. Численное решение сводится к интегрированию системы уравнений с трехдиагональной матрицей коэффициентов. Интегрирование осуществляется с помощью неявной дискретизации по времени, что позволяет получать безусловно устойчивые численные схемы. Конечная система нелинейных алгебраических уравнений допускает эффективное решение с помощью итерационной процедуры Ньютона совместно с обычной прогонкой для систем уравнений с трехдиагональной матрицей. Предложенные алгоритмы являются конечно-элементными и могут быть формально получены с позиций МКЭ. В рамках изложенного подхода наиболее эффективно реализованы возможности МКЭ для построения уравнений, аппроксимирующих задачу в точках разрыва характеристик решения, и построения вычислительных алгоритмов на подвижных сетках, связываемых с особенностями определяемого решения.

В работе подробно рассматривается применение подвижных изотермических расчетных сеток, предоставляющих удобные способы дискретизации в случаях наличия фазовых превращений, неограниченной протяженности области исследования, при интегрировании нелинейных краевых задач, решения которых характеризуются особенностями типа "бегущей волны". Предложен путь обобщения метода изотерм на задачи с немонотонным характером решения. Выполнена аппробация предложенных алгоритмов на большом числе примеров, включающих автомодельные постановки задачи Стефана для плоской и осевой симметрии соответственно при граничных условиях 1-го и И-го рода; замерзание воды в замкнутом объеме с учетом развиваемого давления; задачи с конвективным слагаемым на примере одномерного уравнения Бюргерса при различной величине коэффициента при диффузионном члене от 1 до 10~4 (последний случай дает предельное приближение к решению гиперболического уравнения конвекции); задачу о замерзании раствора соли, на которой проиллюстрировано применение разработанного метода к системам уравнений тепломассопереноса при фазовых превращениях. Сравнение с точными решениями показало, что во всех случаях для достижения погрешности, не превышающей 1 %, достаточно от 3 до 11 узлов изотермической сетки.

Изложенный общий подход допускает и другие способы задания расчетных сеток. Для решения многофронтовых задач Стефана с периодическим

граничным режимом, которые имеют важное значение в проблемах геокриологии, применялась сетка с фиксированными узлами, дополняющаяся по мере возникновения фазовых фронтов соответствующими узлами - изотермами. В качестве примера таких задач рассмотрено решение задачи о "тепловом клапане" (R.Gilpin and B.Wong), моделирующей протаиванне-промерзание поверхностного слоя грунта с учетом сезонного изменения температуры и теплового сопротивления покрова на поверхности в рамках многофронтовой задачи Стефана с граничным условием в форме

-XdTjd х = h(Ja-T), где h = hm + Ahsmait + th \ Ta = Tm +A7sin<u(f + tT) - гармонические функции времени, описывающие суммарный коэффициент теплопередачи к поверхности грунта и температуру воздуха, (о = 2л / tp, tp - величина периода (один год), tb и tT - фазовые смещения, индекс m указывает на среднегодовое значение величины. На рис. 3 представлены два различных результата решения задачи в зависимости от величины среднегодовой температуры воздуха.

1.5

V

0 1гз-ч -^лет 0 < г 3 "^лет

Рис. 3. Протаивание-промерзание при различной среднегодовой температуре воздуха

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

В качестве другого применения разработанного метода к решению многофронтовых задач Стефана отметим численное моделирование эволюции криолитозоны на арктических шельфах (Ю.А.Сигунов, А.И.Фартышев).

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

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

= (9)

илдх (Р.)„ дх Гдт\цъ дг) м,

.дс Р„ дх~дс 1 д f дс\

dt рг дх дг гдг\ дг)

где Р - давление, m - общая пористость, х ■ насыщенность пор льдом, /3(х) и к( х) - эффективная сжимаемость и проницаемость пород, ц - вязкость, рш и рл - плотности воды и льда, vB - скорость фильтрации жидкости, Т -температура, (Ср) = (Ср)(т,х) и Х = Л(т,х) • средневзвешенные по объему пористой среды теплоемкость и коэффициент теплопроводности, L -удельная теплота фазового перехода вода-лед, с - концентрация растворенных солей, D - коэффициент диффузии.

Уравнения (9) - (11) справедливы для областей с непрерывным распределением льдонасьпценности, включая талую зону и область фазовых превращений воды и льда, протекающих в спектре отрицательных температур, не содержащую поверхностей разрыва льдонасьпценности. В последней распределение температуры определяется величиной концентрации солей в не-замерзшем растворе соотношением

T(r,t) = -kc(r,t) при x(r,t) > 0.

Действие механизмов переноса растворенной соли приводит к возникновению скачков льдонасыщенности в зоне фазовых переходов, что требует дополнения основных уравнений (9) - (11) соответствующими условиями на поверхностях разрыва.

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

Исследуемая математическая модель представляет интерес для анализа целого ряда явлений, связанных с фазовыми превращениями в пористых средах, содержащих растворенные примеси, таких как образование ледяных линз в засоленых грунтах или термоэрозионные процессы на морских побережьях, сложенных засоленными частично-мерзлыми породами. Ранее аналогичная система уравнений была исследована в автомодельной постановке в пренебрежении конвективной составляющей потоков тепла и растворенной соли (А.М.Максимов, Г.Г.Цыпкин). В нашей работе главное внимание уделено оценке роли конвективного переноса и методике численного интегрирования систем уравнений тепломассопереноса с подвижными границами скачкообразного изменения льдонасыщенности, возникновение и количественные параметры которых, в отличие от классической задачи Стефана, не заданы в явном виде постановкой задачи. Для прослеживания формирования разрывов используются подвижные сетки, связанные со значениями льдонасыщенности. Возникновение фазовых границ с ненулевым значением скачка \х] фиксируется при слиянии узлов с различной величиной х■ В дальнейшем при слиянии с образовавшейся границей последующих узлов сетки величина скачка на границе корректируется.

Некоторые особенности исследуемого процесса иллюстрируются ниже результатами численных экспериментов. На рис. 4 показано формирование разрывов в распределении льдистости в зоне фазовых переходов в случае инфильтрации жидкости в пласт. Развитый процесс характеризуется в итоге двумя скачками льдистости, представляющими границы частичных фазовых переходов. При этом сохраняется и область непрерывных фазовых превращений, однако их интенсивность в этой области невелика. На рис. 5 показаны распределения безразмерной концентрации соли (сг = с / с, ) вдоль радиальной координаты на различные моменты времени. В начальный период поле концентрации характеризуется наличием "бегущей" волны, продвигаю-

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

ов

Рис. 4

1в

Рис. 5

Аналогичные результаты получены в случае противоположного тока норовой влаги в пласте по направлению к скважине (при отрицательной депрессии па пласт). Перенос соли при течении жидкости из удаленных от скважины пород подавляет процесс разбавления норового раствора в результате таяния льда. Распределение концентрации соли в области охваченной воздействием с течением времени устанавливается близкой к ее значению в норовом растворе в исходном состоянии. Следствием этого является формирование чисто фронтового режима фазовых превращений с исчезновением зоны непрерывных фазовых превращений.

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

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

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

О(0 = {(*.>): gl(y,t)ZxZg2(y,t) , п^У^Уг} в некоторой системе координат, так что положения фазовых границ описываются явными функциями одной из координат и времени х = г„ (у, 0, расчетная сетка определяется семействами линий

х=ги(у,0, и = О, I.....N.

^(У.О = г0(у,О <...< г„(у,г) <...< г^(у.г) = я2(у,г) и У = Ук, к = 0,1-..,К, п =У0 <-<Ук < — УК =Гг-пересекающихся в узлах (л, к.) с координатами (хп,у/с). При этом допускается, что положение линии х = гп (у, О либо задается явно, либо определяется распределением температуры Г = Гп(у,г), так чтобы 7Хг„(>\0.>\ О = Тп(у,г). Во всех случаях считается, что одна из узловых величин (температура Тпк или положение узла гпк на линии у=ук) известна, а другая подлежит определению. Значения ук (к = 0,являются фиксированными. Выполнение оговоренных допущений может быть достигнуто путем подходящего выбора координатной системы. В частности во многих случаях такую возможность дает использование полярной системы координат.

Для внутреннего узла (и, к) его окрестность £)(и,4) представлена объединением элементарных подобластей

¿пЛ) = уди_.14_( 1 где опк = {(Л,У): 2 х < гп+1(у,0, у, < у < ук+1}.-

>,/=0.1

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

(/-/ЧУ + Ю„* = (12)

Ук

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

Суммированием соотношений (12) по всем подобластям окрестности узла с учетом условий баланса тепла на их общих границах (в которых учитывается энергия фазовых превращений, если граница совпадает с поверхностью фазового перехода) выводится тождество, выражающее общий энергетический баланс по окрестности

£(/ + (-1)-,+1 F + У + (-1)' + Рп,к-, + Ря,к = 0 , (13)

7,1=0,1

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

Уравнения относительно узловых неизвестных получаются последовательным раскрытием интегральных слагаемых в (13) по обеим координатам с аппроксимацией температуры линейными приближениями в пределах интервалов г„ (у, г) < х < гп+1 (у, 0 вдоль каждого фиксированного значения координаты у. Для производных по времени используется неявная аппроксимация. В результате получена система девятиточечных алгебраических уравнений с блочно-трехдиагональной матрицей относительно узловых неизвестных для криволинейной подвижной расчетной сетки, в которой пространственное расположение узлов может задаваться как явным образом, так и определяться заданными величинами температур в узлах. Разработан итерационный алгоритм ее решения, использующий частичные итерации по методу Ньютона и матричную прогонку на каждом итерационном шаге.

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

поля вокруг газодобывающих скважин в линейном кусте после 10 лет эксплуатации при расстоянии между скважинами 20 м. С учетом симметрии задача решалась в бесконечной полуполосе при использовании полярных координат на обобщенной температурной сетке размером 8x9 узлов. Удвоение размеров сетки дало полное сопоставление с приведенным результатом.

-30

2.0

г,м Рис. 6

ю

12

14

¡6

18

20

Самостоятельный интерес представляет выполненное моделирование

протаиваяия пластовых льдов, коитакти-

_í_гч з \ $ б 7 в з рующих с малольдистыми породами в

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

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

1. Разработана модификация метода интегрального баланса для приближенного решения одномерных задач Стефана для внутренней и внешней

геометрии задачи и основных типов симметрии (плоской, осевой и сферической) при нелинейных граничных условиях общего вида на возмущающей поверхности. Обоснован способ получения двухсторонних оценок неизвестного закона движения фазовой границы и на его основе оценена погрешность метода, определяемая величиной безразмерных параметров задачи. Эффективность предложенного метода подтверждена тестовыми расчетами, включающими случаи немонотонного движения фронта фазовых переходов.

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

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

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

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

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

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

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

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

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

Основное содержание диссертации опубликовано в следующих печатных работах:

1. Медведский Р.И., Сигунов Ю.А. Численное моделирование повторяющихся тепловых воздействий на толщу мерзлых пород // Изв. СО АН СССР. Сер. техн. наук. 1988. № 21. вып. 6. С. 72-79.

2. Медведский Р.И., Си1унов Ю.А. Метод решения внутренней двухфазной задачи Стефана с нелинейным граничным условием // Теплофизика высоких температур. 1990. Т. 28, № 2. С. 291-300.

3. Медведский Р.И., Сигунов Ю.А. Изменение агрегатного состояния вещества, контактирующего с ограниченным телом // Инж.-физ. журн. 1989. Т. 56, № 1. С. 145.

4. Медведский Р.И., Сигунов Ю.А. О решении одномерных нелинейных задач теплопроводности на изотермической сетке // Журн. вычисл. математики и математической физики. 1989. Т. 29, № 11. С. 1742-1746.

5. Медведский Р.И., Сигунов ЮЛ. Метод численного решения одномерных многофронтовых задач Стефана // Инж.-физ. журн. 1990. Т. 58, ^ 4. С. 681-689.

6. Медведский Р.И., Сигунов ЮЛ. Моделирование воздействия мерзлых пород на нефтегазовые скважины в рамках контактных задач Стефана // Инж.-физ. журн. 1990. Т. 59, № 4. С. 698.

7. Медведский Р.И., Скляр Ю.Г., Си1унов ЮЛ. Влияние ледяных пробок встволе скважин на смятие колонн // Нефтяное хозяйство. 1991. № 11. С. 17-19.

8. Медведский Р.И., Сигунов ЮЛ. Плавление норового льда при твердении тампонажного раствора в заколонном пространстве скважины // Совершенствование технологических способов строительства глубоких разведочных скважип в Западной Сибири. Тр. ЗапСибНИГНИ. Тюмень. 1988. С. 19-29.

9. Медведский Р.И., Сигунов ЮЛ. Обеспечение условий формирования качественного тампонажного камня за кондуктором в интервале морозных глин // Совершенствование технологических способов строительства глубоких разведочных скважин в Западной Сибири. Тр. ЗапСибНИГНИ. Тюмень. 1988. С. 13-19.

10. Сигунов ЮЛ., Фартышев А.И. Исследование эволюции криолито-зоны арктического шельфа методами математического моделирования // Геология и геофизика. 1991. № 8. С. 24-31.