П. Плотник

Явная разностная схема для задачи о замораживании грунта

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

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

            Постановка задачи

    Функции, характеризующие тепловое состояние, и теплофизические свойства грунтов

    Тепловое состояние грунта может быть характеризовано двумя выражающимися друг через друга функциями точки х и времени t – энтальпией (теплосодержанием) единицы объёма грунта h(х, t) и его температурой Т(х, t). Энтальпия h = h(Т) является монотонно возрастающей функцией температуры, имеющей скачок в точке Т = Т* замерзания грунтовой воды.

/*Величина Т* зависит от от степени солоноватости воды и может быть существенно ниже 0 °С*/.

    Обозначим через h_mel и h_fr предельные значения энтальпии при стремлении температуры к Т* со стороны более высоких и более низких температур соответственно. Величина скачка энтальпии в точке замерзания  

        Р = h_mel – h_fr

равна объёмной теплоте замерзания грунта.

   При температурах, отличных от Т*, функция h(Т) дифференцируема, а производная 

       С = dh/dT                                                            (1) 

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

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

        h_fr = 0,

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

   При сделанных допущениях о постоянстве теплоёмкости в пределах каждой из зон функция h(Т) кусочно линейна и определяется равенствами:

  h(Т) = С_fr (Т – Т*)  при T < T*,                               (2)

  h(Т) = Р + С_mel (Т – Т*)  при T > T*,                      (3)

а при Т = Т* величина h может принимать любые значения из интервала 0 ≤ h ≤ Р.

    Функция Т = Т(h), обратная по отношению к функции h = h(Т), является кусочно линейной, причём, в отличие от h(Т), непрерывной, монотонно неубывающей и определённой для всех значений h:

  Т = Т*+ h/С_fr  при h < 0,                                          (4)

  Т = Т* при 0 ≤ h ≤ Р,                                              (5)

  Т = Т*+ (h – Р)/С_mel при h > Р.                              (6)

    Заметим, что свойство однозначной определённости и непрерывности функции Т(h) сохраняются и без предположения о кусочно линейном характере функции h(Т).  

    В соответствии с законом теплопроводности Фурье для вектора плотности теплового потока q = q(х, t) обычно постулируется соотношение:

   q = – k grad T,                                                         (7)

где k – теплопроводность, зависящая, вообще говоря, от температуры.

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

   Полная система уравнений

   Изменение энтальпии грунта со временем определяется дифференциальным уравнением теплового баланса: 

      ∂h/∂t  = – div q.                                                    (8)

   Равенства (7) и (8) вместе с формулами (4)-(6), выражающими температуру через энтальпию, дают полную систему уравнений относительно неизвестных функций h(х, t), Т(х, t)  и q(х, t).

   В классической теории теплопроводности из равенств (7) и (8) исключают плотность теплового потока q, а затем, считая функцию h = h(Т) дифференцируемой и используя определение теплоёмкости (1), получают одно уравнение относительно температуры – уравнение теплопроводности Фурье:

    С ∂Т/∂t = div(k grad T).

   В задачах о промерзании, когда зависимость энтальпии от температуры имеет разрыв, целесообразно не переходить к одному уравнению, а использовать систему (4)-(8) в её исходной форме. 

   Начальные и граничные условия

   Для выделения из бесконечного множества решений системы (4)-(8) требуемого единственного решения, нужно задать дополнительно начальные и граничные условия.

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

   Наличие частной производной по времени от энтальпии ∂h/∂t в уравнении (8) требует задания в начальный момент времени t = t₀ во всей области D начальных значений энтальпии:

      h(х, t₀) = h₀(х) при х є D.                                    (9)

   В качестве граничных условий на поверхности S в данной задаче выступают заданные значения температуры, вообще говоря, зависящие от времени:

      Т(х, t) = Т_S(t) при х є S.                                         (10)

  Интегральный тепловой баланс  

  Проинтегрировав сначала обе части уравнения (8) по области D и воспользовавшись теоремой Гаусса, а затем проинтегрировав результат по времени от t₀ до t, получаем соотношение интегрального теплового баланса:

      H(t) – H(t₀) = U(t),                                                  (11)

где

     H(t) =  ∫ h(х, t) dV, – суммарная величина энтальпии в охлаждаемой зоне в момент времени t,  

                     D                       

                                 

    Q(t) =  ∫ q(х, t)∙n dS, – тепловой поток через поверхность охлаждаемой зоны в момент времени t,            

                    S

 

                    t

    U(t) = ∫ Q(t') dt', – суммарное количество тепла, полученное грунтом в области D за время от t₀ до t.                                        

                   t₀     

    n – единичный вектор внутренней нормали к поверхности S.

    Метод решения задачи

  Общая схема предлагаемого метода

А. В области D, включая граничную поверхность S, строится пространственная сетка, в узлах которой будут вычисляться неизвестные функции.

В. В соответствии с начальным условием (9) определяются значения энтальпии h во всех внутренних (не совпадающих с граничными) узлах сетки.

С. С помощью численного интегрирования энтальпии h по области D вычисляется начальное значение H(t₀) суммарной энтальпии охлаждаемой зоны.

D. Полагается равным нулю начальное значение суммарного количества полученного тепла U(t₀).

Е. Выбирается последовательность моментов времени t_j  = t₀ + j Δt, где Δt – заданный постоянный шаг, в которые будут вычисляться значения неизвестных функций в узлах сетки.  

F. Для каждого момента времени t_j, начиная с j = 0, производятся следующие вычисления:

1.      По заданному значению энтальпии h определяются по формулам (4)-(6) значения температуры Т во внутренних узлах сетки.

2.      Вычисляются значения температуры в граничных узлах сетки в соответствии с граничными условиями (10).

3.      С помощью разностной аппроксимации операции градиента в узлах сетки определяются значения компонент вектора плотности теплового потока q в соответствии с соотношенем (7).

4.      С помощью численного интегрирования по поверхности S определяется значение теплового потока Q(t_j).

5.      Вычисляется значение суммарного количества тепла U(t_j), полученного грунтом к моменту времени t_j, прибавкой величины Q(t_j) ∙ Δt  к значению на предыдущем шаге.

6.      С помощью разностной аппроксимации операции дивергенция  во внутренних узлах сетки определяются значения дивергенции вектора q.

7.      С помощью разностной аппроксимации дифференцирования по времени в уравнении (8) определяются значения энтальпии h во всех внутренних узлах сетки для следующего момента времени t_j+1.

8.      С помощью численного интегрирования по области D вычисляется значение H(t_j+1) суммарной энтальпии охлаждаемой зоны в момент времени t_j+1.

9.      Проверяется выполнение теплового баланса (11). /*Нарушение баланса означает наличие ошибки в программе*/

10.  По изменению знака энтальпии определяется положение фронта промерзания грунта на следующем шаге.

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

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

   Плоский осесимметричный характер полей

   Рассмотрим замораживание влажного грунта с помощью холодильной колонки – вертикальной трубы с внешним радиусом r₀, внутри которой циркулирует хладагент, имеющий в установившемся режиме температуру Т `.

   Ставится задача найти изменение температурного поля и закон продвижения фронта замерзания грунта в предположении о том, что отсутствуют другие тепловые стоки или источники, а также любые неоднородности среды, способные исказить создаваемое под влиянием колонки осесимметричное температурное поле. Если при этом ограничиться рассмотрением температурных полей на расстояниях от оси колонки, значительно меньших толщины пластов грунта с отличающимися свойствами, то в пределах каждого пласта в цилиндрической системе координат (r, φ, z), ось Oz которой совпадает с осью колонки, можно считать температуру и энтальпию грунта зависящими только от координаты r и от времени t:

   h =h(r, t),  Т = Т(r, t).

   У вектора плотности теплового потока q в этих условиях будет иметься только одна отличная от нуля компонента q, направленная по радиусу, причём

   q = q(r, t).

   Граничные условия 

   Роль области D в рассматриваемом случае играет цилиндрическое кольцо, высота которого L равна толщине рассматриваемого геологического пласта с однородными свойствами, внутренний радиус равен r₀, а величина внешнего радиуса R выбирается настолько большой, что на рассматриваемом интервале времени охлаждение не оказывает заметного влияния на температуру Т(R, t).

   Поверхность S состоит из четырёх частей:

   S = S₁ + S₂ + S₃ + S₄,

где

 S₁ – круговой цилиндр с радиусом r₀ и высотой L, совпадающий с участком поверхности охлаждающей колонки, лежащим в пределах рассматриваемого пласта;

 S₂ – круговой цилиндр с радиусом R и высотой L, на поверхности которого можно считать, что температура грунта неизменна в течение всего процесса замораживания;

 S₃ и S₄ – плоские торцевые кольца с радиусами r₀ и R, на которых не требуется задание граничных условий, т.к. по предположению неизвестные функции не зависят от координаты z.

   Граничное условие на S₁

    Если обозначить через Т₀ температуру грунта перед началом процесса замораживания, то изменение во времени температуры поверхности колонки может быть вычислено по формуле:

    T(r₀, t) = Т₀ exp[(t₀ – t) /τ] + Т `{1 – exp[(t₀ – t)/τ]},    (12) 

где через τ обозначено характерное время выхода холодильной системы на стационарный режим работы.

   Граничное условие на S₂

    T(R, t) = Т₀.                                                                  (13)

    Разностная сетка

    Разобъём отрезок r₀ ≤ r ≤ R на N равных частей длиной

             Δr = (R – r₀) / N

 и обозначим через r_i узлы сетки (границы разбиения):

           r_i = r₀ + i Δr   (i = 0, 1, ..., N).                               (14)

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

          t_j = t₀ + j Δt   (j = 0, 1 , ...),

где Δt – заданный временной шаг.   

   Для значений в точке сетки r_i неизвестных функций, рассматриваемых в один и тот же момент времени t_j, будут использоваться обозначения:

     h_i =h(r_i, t_j), Т_i = Т(r_i, t_j),  q_i = q(r_i, t_j).

   Значение энтальпии в той же точке сетки r_i, но на следующем временном шаге (при t = t_j+1) будет отмечаться значком “^”:

       h_i^=h(r_i, t_j+1).

   Для сеточного значения дивергенции вектора q будет использоваться обозначение:

       d_i = div q при r = r_i при t = t_j.

    Аппроксимация градиента температуры и плотности теплового потока

    У вектора градиента температуры в рассматриваемом случае отлична от нуля только r-компонента, равная ∂Т/∂r. Для определения этой компоненты в узле r_i будет использоваться аппроксимация правой односторонней разностью:

   ∂Т/∂r = (Т_{i +1} – Т_i) / Δr.                                                 (15)

    Из (15) и (7) следует формула для вычисления в узлах сетки i = 0, 1, ..., N– 1 разностной аппроксимации единственной отличной от нуля компоненты вектора плотности теплового потока:

    q_i =  k_i (Т_i – Т_{i +1}) / Δr.                                                  (16)   

где

   k_i = k_fr  если Т_i ≤ Т*,

   k_i = k_mel в противном случае.

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

    Аппроксимация дивергенции вектора плотности теплового потока

    В формуле для дивергенции вектора q в цилиндрических воординатах фигурирует только отличная от нуля компонента q = q(r, t):  

    div q = r ^–1 ∂(r q) / ∂r.

    Используя в этом случае для частной производной по r аппроксимацию левой односторонней разностью, имеем:

    d_i = r_i ^–1 (r_i q_i. – r_{i –1} q_{i –1}) / Δr.                                     (17)   

    Формула (17) позволяет определить дивергенцию плотности теплового потока в узлах сетки i = 1, ..., N– 1.

    Аппроксимация частной производной по времени от энтальпии

   ∂h/∂t (r_i, t_j) = (h_i^ – h_i) / Δt.                                          (18)

    Определение энтальпии на новом временном слое

      С учётом (17) и (18) разностный аналог дифференциального уравнения локального теплового баланса (8) имеет вид:

      (h_i^ – h_i) / Δt = r_i ^–1 (r_{i –1} q_{i –1} – r_i q_i) / Δr.                  (19)   

      Из (19) непосредственно следует формула, по которой могут быть рассчитаны значения энтальпии в момент времени t = t_j+1 в тех же точках, в которых была определена дивергенция вектора q, т.е при i = 1, ..., N– 1.

       h_i^ = h_i + Δt (r_i ^–1r_{i –1} q_{i –1} – q_i) / Δr.                          (20)

    Тепловой баланс на разностной сетке

      Умножая обе части (19) на произведение 2π L r_i Δr Δt и суммируя полученный результат по i от 1 до N – 1, получаем соотношение

       H(t_j+1) – H(t_j) = Q₀ – Q_N–1,                                         (21) 

   где

                                     N –1

     H(t) = 2π L Δr Σ r_i h_i(t)                                               (22)                                           

                                     i = 1

       – суммарная энтальпия охлаждаемой зоны в момент времени t;

     Q₀ = 2π L r₀ q₀ Δt                                                         (23)                                                 

     – теплоприток от холодильной колонки на интервале времени от t = t_j до t = t_j+1;   

     Q_N–1 = 2π L r _{N – 1} q _N–1 Δt                                             (24)

       – теплоприток от неохлаждённой зоны на интервале времени от t = t_j до t = t_j+1.     

/*На самом деле оба «теплопритока» Q₀ и Q_{N – 1} отрицательны и являются теплоотводами*/           

     Просуммируем теперь обе части равенства (21) по временным шагам от t = t₀ до t = t_j+1. В результате получаем соотношение интегрального теплового баланса, являющееся разностным аналогом равенства (11):

       H(t_j+1) – H(t₀) = U₀ (t_j+1) – U_N–1 (t_j+1),                         (25) 

где U₀ (t) и U_N–1 (t) – суммарные теплопритоки за всё время, прошедшее от начала процесса до момента времени t, соответственно от холодильной колонки и от неохлаждённой зоны.

     Величины U₀ (t) и U_N–1 (t) могут быть вычислены по рекуррентным формулам:

      U₀ (t _j+1) = U₀ (t _j) + Q₀,                                                (26)

      U _N–1 (t _j+1) = U _N–1 (t _j) + Q _N–1.                                     (27)

    Алгоритм расчёта

 Предварительный расчёт до начала временного цикла

1.      Объявляются массивы вещественных переменных r_i, h_i, Т_i, q_i  (i = 0, 1, ..., N).

2.      По формуле (14) определяются элементы массива r_i для i = 0, 1, ..., N.

3.      По формуле (3) и заданному начальному значению температуры Т = Т₀ вычисляются элементы массива h_i во внутренних точках (для i =1, ..., N – 1).   

4.      Вычисляется начальное значение H(t₀) суммарной энтальпии охлаждаемой зоны по формуле (22).

5.      Полагаются равными нулю начальные значения суммарных теплопритоков U₀ и U_N–1. 

 Расчёт временного цикла

 Для каждого момента времени t_j, начиная с j = 0, производятся следующие вычисления:

1.      Для i =1, ..., N – 1 по известным значениям энтальпии h_i определяются по формулам (4)-(6) значения температуры Т_i.

2.      По формулам (12) и (13) вычисляются значения температуры Т_i соответственно для i =0 и  i = N.

3.      По формуле (16) вычисляются значения q_i для i = 0, 1, ..., N – 1.

4.      По формулам (23) и (24) определяются величины теплопритоков Q₀ и Q_N–1.

5.      По формулам (26) и (27) вычисляются значение суммарного количества тепла U₀ (t _j+1) и U _N–1 (t _j+1), полученного грунтом к моменту времени t_j+1 соответственно от холодильной колонки и от неохлаждённой зоны.

6.      По формуле (20) определяются для i =1, ..., N – 1 значения энтальпии h_i в момент времени t _j+1.

7.      По формуле (22) вычисляется значение H(t_j+1) суммарной энтальпии охлаждаемой зоны в момент времени t_j+1.

8.      Проверяется выполнение условия теплового баланса (25). /*Баланс должен выполняться с большой точностью (обычно невязка оказывается машинным нулём). Нарушение баланса означает наличие ошибки в программе*/

9.      Последовательным перебором элементов массива h_i определяется значение радиуса s фронта замерзания в момент времени t_j+1, как максимально значение r_i, для которого выполняется условие h_i ≤ 0.

10.  Принимается решение о продолжении или окончании временного цикла по условиию достижения требуемого положения s* фронта промерзания или по истечению контрольного срока времени.

     Исходные данные

С_mel, Дж/(м³·К) – теплоёмкость единицы объёма талого грунта;

С_fr, Дж/(м³·К) – теплоёмкость единицы объёма мёрзлого грунта;

k_mel, Вт/(м · К) – коэффициент теплопроводности талого грунта;

k_fr, Вт/(м · К) – коэффициент теплопроводности мёрзлого грунта;

L, м – толщина рассматриваемого геологического пласта с однородыми теплофизическими свойствами;

Р, Дж/м³ – скрытая теплота замерзания единицы объёма грунта;

R, м – расстояние от оси охлаждающей колонки, на котором в течение всего процесса замораживания существенного понижения температуры не происходит; /*Подбирается в ходе пробных расчётов из условия, чтобы во все моменты времени модуль теплопритока со стороны неохлаждаемой зоны Q_N–1 составлял не более 1%  модуля теплопритока от холодильной колонки Q₀.*/  

r₀, м – радиус холодильной колонки;

s*, м – требуемое конечное положение фронта замерзания грунта;

Т*, °С – температура замерзания грунта;

Т `, °С – температура стенки холодильной колонки в установившемся режиме;

Т₀, °С – температура грунта перед началом процесса замораживания;

Δr, м – пространственный шаг разностной сетки;

Δt, с – временной шаг;

/*Величины шагов Δr и Δt должны удовлетворять условиям устойчивости разностной схемы:

  k_mel Δt /(С_mel Δr ²) ≤ ½ ; k_fr Δt /(С_fr Δr ²) ≤ ½.

   Обычно эти условия выполняются, если принять:

 Δr = 0.1 м;

 Δt = 3600 с. */

τ, с  – ориентировочное время выхода на режим холодильной установки.

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

Tаблица 1. 

____________________________________________________________

                                  Мокрый песок     Мокрая глина          Вода/лёд

____________________________________________________________

k_mel, Вт/(м · К)                    1.9                         1.4                        0.6       

С_mel, Дж/(м³·К)                 2.8∙10⁶                    3∙10⁶                   4.2∙10⁶

k_fr, Вт/(м · К)                      2.1                         1.6                        2.3

С_fr, Дж/(м³·К)                   2.1∙10⁶                    2.2∙10⁶                1.9∙10⁶

Р, Дж/м³                           6.7∙10⁷                   1.05∙10⁸              3.32∙10⁸

                                                                                                                                                                           Дата последних изменений: 2011-07-13

Главная страница