Численное моделирование процессов тепломассопереноса в испаряющемся слое жидкости

Для решения задачи (7) – (12)  заменим область непрерывного изменения аргументов x и t областью дискретного их изменения. Так, отрезок изменения аргумента x разделим на N одинаковых частей с шагом

h = 1/N. Аналогично поступим с аргументом t. В качестве шага по времени t возьмем τ. Таким образом, получим сетку, образованную пересечением прямых x_n и t_m , n = 0,..,N, m = 0, 1, … Сеточные функции, определенные в узлах сетки (x_n ,t_m) являются дискретными аналогами искомых непрерывных функций [5].

Качество сеточной схемы  в основном определяется выбранным  шаблоном, который задает номера узлов для аппроксимации производных  
Чисто неявной разностной схемой для уравнений (7) – (8) называется  разностная схема, которая использует 4-х точечный шаблон вида:

 

 (x_n-1,τ_m+1)        (x_n,τ_m+1)           (x_n+1,τ_m+1)          

 

 

 

(x_n,τ_m)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.1.1. РАЗНОСТНЫЕ  СООТНОШЕНИЯ ДЛЯ УРАВНЕНИЯ

ТЕПЛОПРОВОДНОСТИ

Построим неявную разностную схему для уравнения теплопроводности (7)  в соответствие с выбранным шаблоном.

Разностные аппроксимации частных  производных первого и второго  порядка имеют вид:

где – значение сеточной функции в узле (x_m,τ_n). Тогда для уравнения теплопроводности получаем:

                                               (14)     

Граничные и начальные условия имеют вид:

                                                           (15)           

Эта схема имеет порядок  аппроксимации  .

Каждое разностное уравнение (14) содержит на каждом новом слое три неизвестных значения:  ,  ,  , которые нельзя сразу определить через известное решение на предыдущем m-том слое. Такие схемы называются неявными. Заметим, что разностная схема (14) состоит из линейных алгебраических уравнений с трехдиагональной матрицей и для ее решения можно применить метод прогонки.

Преобразуем полученную схему, чтобы  можно было применить ее к решению  методом прогонки [5]:

                                  (16)

Получили систему из N+1 линейных алгебраических уравнений на (m+1) шаге по времени.

Выпишем матрицу для этой системы  с учетом начальных и граничных  условий:

                             (17)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.1.2. РАЗНОСТНЫЕ  СООТНОШЕНИЯ ДЛЯ УРАВНЕНИЯ ДИФФУЗИИ

 

Построим разностную схему для уравнения диффузии (8).

Разностные аппроксимации частных  производных первого и второго  порядка имеют вид:

,

где – значение сеточной функции в узле (x_m,τ_n). Тогда для уравнения диффузии получаем:

                            (18)

Граничные и начальные условия  принимают вид:

                                                                 (19)

Эта схема имеет порядок аппроксимации  .

Преобразуем полученную схему, чтобы  можно было применить ее к решению  методом прогонки [5]:

                               (20)

Получили систему из N+1 линейных алгебраических уравнений на (m+1) шаге по времени.

Выпишем матрицу для  этой системы с учетом начальных  и граничных условий:

                             (21)

 

2.1.3. МЕТОД ПРОГОНКИ

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

Преобразуем разностную схему к  следующему виду:

                                            (22)

с граничными условиями вида:

                                                                                     (23)

                                                                                   (24)

Будем использовать решение из (13)-(15) в следующей форме:

                                                 (25)

где:

                                                                                     (26)

                                                              (27)

α₁ и  β₁ находим из (25) и краевого условия (23):

Значение y_p , необходимое для начала счета по формулам (25), находим из (25) и краевого условия (24) при :

                                                                                     (28)

Матрица системы имеет вид:

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

 

 

 

 

 

 

2.1.4. АППРОКСИМАЦИЯ УРАВНЕНИЯ ДВИЖЕНИЯ ГРАНИЦЫ.

 

Смещение границы слоя раствора считаем с помощью явного метода Эйлера. Тогда уравнение (9) примет вид:

                                                    (29)

где Q < 0, поскольку направление  поступления тепла в слой противоположно положительному направлению оси Ox.

Начальное условие имеет вид:                                                     (30)

 

2.1.5. МЕТОД ЭЙЛЕРА

Для решения задачи Коши используется явный метод Эйлера . Данный метод  является универсальным методом приближенного решения задачи Коши  y'(x) = f(x,y(x))                                                                                         (31)

            y(x₀)=y₀.                                                                                                    (32)

Под решением такой задачи понимается функция y(x), заданная на отрезке [x₀, x₀+X] (X>0, X – длина отрезка интегрирования), обладающая в каждой точке этого отрезка производной, удовлетворяющая в каждой точке отрезка уравнению (31), а при x = x₀ – дополнительному («начальному») условию (32). Существование и единственность такого решения мы, естественно, предполагаем. Кроме того, чтобы гарантировать существование и приближенного решения, условимся считать, что правая часть F(x, y(x) )  уравнения (31) определена в любой точке (x_* , y_*) полосы, заданной на плоскости переменных x, y неравенствами x₀ ≤ x ≤ x₀+X  (рис. 1)

Зададимся натуральным  N  и разделим отрезок интегрирования [x₀, x₀+X] на N равных частей длины h = X / N  (33)  точками x_i = x₀ + i*h , i=0,1, . . . ,N  (34). Дискретное множество точек (34) называется сеткой на отрезке [x₀, x₀+X], а сами точки x_i – узлами этой сетки. Расстояние между узлами сетки с соседними номерами равно длине (33) частичного отрезка разбиения [x_i, x_i+1], эта величина называется шагом сетки. При неограниченном увеличении N шаг сетки стремится к нулю:  h→0  при N→∞, и сетка становится все более и более густой.

 В общем, цель состоит в  том, чтобы вывести систему  уравнений для приближенного отыскания значений y(x_i) искомого решения y в узлах сетки. Для этого следует сделать замену в дифференциальном уравнении (31), записанном в узле сетки x_i : y' (x_i) = f (x_i , y(x_i))                                         (35)

производной y'(x_i) ее простейшим сеточным приближением

В методе Эйлера приближенные значения у(х_i) » y_i вычисляются последовательно по формулам у_i+h * f(x_i, y_i) (i=0,1,2…).При этом искомая интегральная кривая у=у(х), проходящая через точку М₀(х₀, у₀), заменяется ломаной М₀М₁М₂… с вершинами М_i(x_i, y_i) (i=0,1,2,…); каждое звено М_iM_i+1 этой ломаной, называемой ломаной Эйлера, имеет направление, совпадающее с направлением той интегральной кривой уравнения (31), которая проходит через точку М_i. Если правая часть уравнения (31) в некотором прямоугольнике R{ |x-x₀| £ a, |y-y₀| £ b }удовлетворяет условиям: |f(x, y₁)- f(x, y₂)| £ N|y₁-y₂| (N=const),|df/dx|=|df/dx+f(df/dy)| £ M (M=const), то имеет место следующая оценка погрешности: |y(x_n) - y_n| £ h * M /  2 *N*[ (1+h*N)ⁿ - 1] (36), где у(х_n)-значение точного решения уравнения (31) при х=х_n, а у_n- приближенное значение, полученное на n-ом шаге. Формула (36) имеет в основном теоретическое применение.

Метод Эйлера легко распространяется на системы дифференциальных уравнений  и на дифференциальные уравнения высших порядков. Последние должны быть предварительно приведены к системе дифференциальных уравнений первого порядка [6].

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.2. АЛГОРИТМ  РЕШЕНИЯ ЗАДАЧИ

 

Поставленная задача с учетом разностных схем принимает вид:

,

,

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

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

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

Следуя заданным предположениям, начинаем решение задачи с вычисления поля температуры по формулам (16) методом прогонки. Граничное условие вычисляем по формуле (24), подставляя туда соответствующие значения и полагая, что концентрация еще не изменилась и равна начальным значениям. Затем вычисляем новое положение границы по формуле (29). Переносим границу в новый узел, если требуется. Далее вычисляем новые значения поля концентрации растворителя по формулам (20) методом прогонки. Если после расчета поля температуры, толщины слоя, граница была перенесена в новый узел, то значение температуры берем в соответствии с предыдущим положением границы.

Последовательность действий при  решении задачи:

1. Расчет нового значения температуры.
2. Расчет нового положения границы.
3. Перенос границы в новый узел, если требуется.
4. Расчет концентрации.

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3.  ПРОГРАММНАЯ РЕАЛИЗАЦИЯ ПРЕДЛОЖЕННОГО                  

                                 РЕШЕНИЯ ЗАДАЧИ

 

    3.1. РАЗРАБОТКА ПРОГРАММЫ.

Для реализации алгоритма  численного решения задачи был выбран язык программирования Delphi 7. Были разработаны методы вычисления метода прогонки, метода Эйлера, вычисления нового положения границы. А также реализована заявленная в п. 2.2 последовательность действий для решения задачи. Результаты выводятся на графиках. Вычисление поля температуры реализовано в конструкторе Tprogon.Create, концентрации в конструкторе TprogonK.Create.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

                  3.2 ТЕСТИРОВАНИЕ ПРОГРАММЫ

 

Проверим работу программы  на задаче, решение которой известно.  

Рассмотрим условие, при  котором испарение отсутствует, т.е. Q=0, R=0, α_D=1. Тогда задача (7) – (9) перепишется в виде

                                                                                                       (37)

                                                                                                          (38)

.                                                                                                      (39)

Начальные и граничные  условия для T и k остаются такими же, как (10) и (11), но T_δ и k_δ теперь заданные константы, а граничное условие для k при x=0 равно k(0,t) = k₀.

При [5] .

Находим решение этой задачи:

Из краевых условий  Отсюда:

 из условия задачи. Тогда при

Рис. 1

 

При 

Рис. 2

         3.3. ЧИСЛЕННЫЙ АНАЛИЗ ЗАДАЧИ