ПОФП учебник Кошмаров справленный

6.12. Метод численного решения дифференциальных уравнений

6.12.1. Замкнутая система уравнений приведенной в данном пособии математической модели может быть решена только численным способом. Современный уровень развития вычислительной техники и вычислительной математики позволяет разработать достаточно точные численные методы, в которых реализуются различные конечно-разностные схемы исходной системы уравнений. Используем наиболее апробированный численный метод решения систем дифференциальных уравнений в частных производных - метод «контрольных объемов» [15], используемый в современных программных пакетах расчета на ПЭВМ прикладных задач газодинамики и тепломассообмена PHOENICS [21] и FLUENT [22].

6.12.2. Дифференциальные уравнения в частных производных, приведенные к «стандартному» виду (1.10), решаются конечно-разностным методом "контрольных объемов" по неявной схеме [15]. Для трехмерной задачи дискретный аналог уравнений имеет следующий вид (рис.6. 4.1):

(6.12.1)

где ; ; ; ; ; ; точки T и В расположены аналогично в плоскостях YOZ и XOZ; ; ; ; – значение плотности на предыдущем шаге по времени, кг/м³; - значение функции на предыдущем шаге по времени; Dx, Dy, Dz – шаги вдоль осей OX, OY, OZ соответственно, м; Dt - шаг по времени, с.

6.12.3. Используется функция в следующем виде (схема со степенным законом изменения числа Пекле внутри контрольного объема):

; (6.12.2)

Операторы типа и означают, что при расчете используется большая из величин, разделенных запятой.

б в

Рис. 6.12.1. Шаблоны контрольных объемов: а - для основных узлов конечно-разностной сетки; б - для проекции скорости на ось ОХ; в - для проекции скорости на ось ОУ

6.12.4. Число Пекле Р в выражении (4.2) определяется как отношение F и D, таким образом, P_e=F_e/D_e и т. д. Расходы F и проводимости D определяются следующим образом:

(4.3)

6.12.6. Учитывая, что при пожаре плотность газовой среды помещения существенно изменяется, необходимо использовать уравнение для определения поправки давления в «сжимаемой» форме:

, (6.12.4)

где p_i^’ – поправки давления в соответствующих точках.

6.12.7. Поправочные формулы для проекций скорости имеют следующий вид:

; (6.12.5)

; (6.12.6)

, (6.12.7)

где величины с верхним индексом ^* означают текущие значения величины, а с индексом ¢ – поправочные значения; d_e, d_n, d_t – величины, получаемые из соответствующих уравнений для компонентов скорости.

6.12.8. Расчёты проводятся на неравномерных конечно-разностных сетках (Dx=var, Dy=var, Dz=var) с переменным шагом по времени, определяемым из условия Куранта [15] для всех узлов конечно-разностной сетки:

, (6.12.8)

где w_a – локальная скорость звука, м/с, k_t<1 – число Куранта.

6.12.9. Последовательность решения системы алгебраических уравнений, приведенных к «стандартному» виду (6.12.1), является следующей:

1) задается приближенное поле давления (например, давление во всех контрольных объемах равно атмосферному);

2) решаются алгебраические уравнения (6.12.1), соответствующие уравнениям движения (6.2)-(6.4), для получения полей значений скоростей w_x^*=f₁(x,y,z), w_y^*=f₂(x,y,z)и w_z^*=f₃(x,y,z);

3) решается уравнение (6.12.4) для определения поправки давления p^’=f₄(x,y,z);

4) рассчитывается скорректированное поле давлений p=р^*+р';

5) проводится расчет составляющих скорости по формулам (4.5)- (4.7);

6) находится решение алгебраических уравнений (6.12.1), соответствующих уравнениям неразрывности смеси (6.1) и ее отдельных компонентов (6.6), уравнению энергии (6.5), уравнений теплопроводности (6.11) - (6.14), а также уравнений (6.19) и (6.20) k-e модели турбулентности, для определения полей других физических величин (температуры газовой среды и строительных конструкций, концентрации компонентов смеси газов, коэффициентов тепломассопереноса, теплофизических свойств);

7) представляется скорректированное давление p как новое p^* и расчет продолжается с пункта 2.

Процедура расчета повторяется до тех пор, пока не будет получено сходящееся решение, т.е. отличие параметров газовой среды на соседних итерациях не должно превышать заранее заданного значения (например, ½Т⁽^k⁺¹⁾-T⁽^k⁾½<0,01 К, где k+1, k – номера текущей и предыдущей итераций).

6.12.10. Используется неявная конечно-разностная схема [15]. При этом система линейных алгебраических уравнений (6.12.1) решена методом продольно-поперечной прогонки [15]. Сначала производится решение для одного из выбранных направлений (например, вдоль оси ОХ) с помощью стандартного метода исключения Гаусса. Каждое уравнение системы приводится к следующему виду:

, (6.12.9)

где i=1, 2,..., N_x – номер точки конечно-разностной схемы вдоль оси ОХ.

Таким образом, физическая величина Ф_i связана с соседними значениями этой величины Ф_i-₁ и Ф_i+₁. Далее предполагается, что существует следующая зависимость:

, (6.12.10)

где a_i и b_i – прогоночные коэффициенты.

Подставляя формулу (6.12.10) в (6.12.9), получаем зависимость прогоночных коэффициентов a_i и b_i от a_i-₁ и b_i-₁:

; (6.12.11)

. (6.12.12)

Значения a₁ и b₁ определяются исходя из граничных условий:

; . (6.12.13)

Например, в случае граничных условий второго рода при прогонке вдоль оси ОХ для уравнения энергии (6.5) можно записать:

, (6.12.14)

где – градиент температуры по нормали к поверхности ограждающих конструкций на ее границе; Т₁, Т₂ – температуры в центрах граничного и соседнего контрольных объемов; Dx – шаг вдоль оси ОХ.

Тогда, в соответствии с уравнениями (4.10) и (4.14) значения a₁ и b₁ равны:

; . (6.12.15)

После этого, все прогоночные коэффициенты рассчитываются по формулам (6.12.11) и (6.12.12).

Далее исходя из граничных условий на другой границе (i=N_x) определяется величина Ф_Nx и по формуле (6.12.10) все остальные значения Ф_i. Например, при граничных условиях второго рода в случае прогонки вдоль оси ОХ (уравнение энергии (6.5)):

; , (6.12.16)

где Т_Nx, Т_Nx_-1 – температуры в центре граничного к поверхности и соседнего контрольных объемов.

Затем выполняется аналогичная прогонка вдоль двух других осей ОУ и ОZ с использованием в качестве значений Ф_i последних в итерационном смысле величин. Описанный алгоритм повторяется вдоль всех осей до тех пор, пока значения функции Ф_iна соседних итерацияхво всех узлах конечно-разностной сетки будут отличаться друг от друга на заранее заданную величину.

6.12.11. Сложная геометрия конструкций и громоздких предметов в помещении задается блокированием в области твердых поверхностей контрольных объемов конечно-разностной сетки, используемой при численном решении полевой модели на персональном компьютере, с помощью задания завышенной величины вязкости смеси газов: m=10¹⁴ кг/(м×с).

Содержание