3.2 Расчет турбулентного течения в трубах
Задание – рассчитать инерционные силы и силы вязкости, действующие на поверхность канала течения при турбулентном движении промывочной жидкости по внутритрубному пространству.
Для решения задачи используются следующие исходные данные
Расчетным элементом является канал
течения промывочной жидкости длиной
,
представляющий собой внутреннюю полость
бурильной трубы типа ПН-12710.
Реологическими характеристиками промывочной жидкости, определенными по результатам реометрии, являются:
- верхняя граница средней скорости
сдвига, соответствующая началу
турбулентного течения –
=240
с-1;
- нижняя граница средней скорости сдвига,
соответствующая началу пластического
течения –
=16
с-1;
- пластическая вязкость –
=0,4779
Па*с;
- реологическая модель Carreau, построенная по точкам реологического профиля, лежащим в области ламинарного течения, и пластической вязкости
;
- модель удельного повышения потерь давления при турбулентном течении
.
Плотность промывочной жидкости составляет пж=1120 кг/м3, ее расход при промывке скважины – Q=0,056 м3/с.
Решение
Определяем удельные потери давления,
возникающие при ламинарном течении
бурового раствора внутри труб, используя
методику предыдущей работы. В результате
вычислений получаем: а) удельные потери
давления, равные
=8338,537
Па/м; б) среднюю скорость сдвига
в сечении потока
=310,42
с-1.
Рассчитываем среднюю скорость движения жидкости вдоль оси трубы
м/с. (3.1)
Рассчитываем осредненную эффективную
вязкость ламинарного потока, соответствующую
данной скорости сдвига, подставляя
в уравнение реологической модели
(3.2)
Вычисляем перепад давления на торцах
расчетного элемента при турбулентном
режиме течения, принимая длину расчетного
элемента
=1
м. Для этого используем формулу
, (3.3) где 
.
Получаем
Па.
Создаем модель в системе мультифизического моделирования «FEMLAB».
Выбираем глобальные параметры модели. Для этого, открыв список моделей (Model Navigator), выбираем турбулентный поток несжимаемой жидкости в модуле "Chemical engineering module", как показано на рисунке 3.15.
Рисунок 3.15 – Выбор модели с использованием навигатора
Как показано на рисунке 3.16, в режиме "Draw Mode" задаем геометрию расчетного элемента в виде прямоугольника, соответствующего половине продольного сечения канала течения жидкости.
Рисунок 3.16 – Задание геометрии модели
В качестве ширины прямоугольника
указываем внутренний радиус трубы.
Длина расчетного элемента при моделировании
турбулентного потока должна обеспечивать
воспроизведение пульсаций как в малом,
так и большом масштабе времени, причем
желательно перекрыть интервал нескольких
крупномасштабных периодов. Помимо этого
поток, заданный на входе в элемент
средней скоростью, должен установиться
в соответствии с турбулентным профилем,
что происходит на определенном расстоянии
от входа. В этой связи длину элемента
имеет смысл выбирать настолько большой,
насколько это позволяют возможности
ПЭВМ. В нашем случае расчетные возможности
программы ограничены имеющейся
оперативной памятью, размер которой
допускает число элементов сетки на
уровне 8000-10000. Этому соответствует дина
элемента около одного метра. Для упрощения
интерпретации результатов расчета
программы длину целесообразно выбрать
целым числом, поэтому останавливаемся
на ее величине
=1,0
м.
Задаем граничные условия, как это показано на рисунке 3.17, перейдя в режим "Boundary Mode", для чего выбираем команду [Boundary settings...] на вкладке [Physics] главной панели инструментов «FEMLAB».
Рисунок 3.17 – Задание граничных условий
На входе в элемент (граница №2) задаем
значения составляющих скорости жидкости
по осям координат. При этом в качестве
составляющей скорости вдоль оси течения
(z) используем значение
,
вычисленное ранее.
Задаем свойства сплошной среды (жидкости) между границами области расчета в соответствии с рисунком 3.18. Для этого переходим в режиме "Subdomain Mode", используя команду [Subdomain settings...] на вкладке [Physics] панели инструментов «FEMLAB».
Рисунок 3.18 – Задание физических свойств сплошной среды
В поле {Kinematic viscosity} на вкладке [Physics] окна «Subdomain settings – k-epsilon turbulence model (ke)» задаем отношение вычисленной ранее осредненной эффективной вязкости жидкости к ее плотности, которая также указывается независимо в поле {Density} этой вкладки. Затем щелчком левой кнопки мыши по кнопке <Artificial Stabilization>, расположенной на той же вкладке, активируется соответствующее окно, в котором необходимо задействовать опции {Isotropic diffusion} и {Streamline diffusion}, и установить значения параметров их настройки (Tuning parameter), как это показано на рисунке 3.18.
Разбиваем область расчета на конечные элементы, для чего переходим в режим “Mesh Mode”, используя команду [Mesh Parameters…] на вкладке [Mesh] главной инструментальной панели «FEMLAB». При этом открывается окно задания параметров сетки (Mesh Parameters), на вкладке [Global] которого указывается максимальный размер элемента в поле {Maximum element size:}, как показано на рисунке 3.19. При этом исходим из необходимости обеспечения количества конечных элементов вдоль границы области расчета, соответствующей поперечному сечению потока, не менее 10.
Рисунок 3.19 – Задание параметров сетки и результат ее построения
В результате получается сетка, количество элементов которой, согласно информационному табло внизу рабочего окна системы, показанного на рисунке 3.19, составляет 9066.
Устанавливаем параметры решателя, как это показано на рисунке 3.20, используя кнопку <Solver parameters> на панели инструментов «FEMLAB» для открытия соответствующего окна.
На вкладке [General] данного окна в поле {Name of parameter:} указываем имя переменной, соответствующей давлению на выходе из элемента (pout), которое должно совпадать с именем аналогичной переменной, используемым при задании граничных условий (рисунок 3.17).
Рисунок 3.20 – Установка параметров решателя
В поле {List of
parameter values}
задаем два крайних значения переменной
pout, для которых
будет проведен расчет. Указанные значения
получаем следующим образом. Верхний
предел (
)
– путем увеличения
на 15 % и округления результата в большую
сторону. Нижний предел (
)
– путем уменьшения
на 5 % и округления результата в меньшую
сторону.
Решаем задачу в следующей последовательности.
Выполняем расчет при крайних значениях переменной pout, используя кнопку <Solve problem> на панели инструментов «FEMLAB» или посредством инструмента "Solver manager". В результате получаем набор распределений фазовых переменных по области решения. Результаты расчета для фазовой переменной – поле скоростей (Velocity Field) показаны на рисунке 3.21.
Рисунок 3.21 – Решения задачи для фазовой переменной "поле скоростей"
при давлении на входе в элемент 10000 и 8000 Па
Определяем значения давления на входе в элемент, используя интегратор, активируемый посредством меню, как изображено на рисунке 3.21. Значения интегралов отображаются на информационном табло внизу рабочего окна «FEMLAB» и обозначены на рисунке выделением. Для получения значения входного давления необходимо разделить результат интегрирования на площадь поперечного сечения канала
.
(3.4)
Таким образом, получаем:
- для
=10000
Па
Па;
- для
=8000
Па
Па.
Вычисляем перепад давления на элементе
.
(3.5)
Получаем:
- для
=10000
Па
Па;
- для
=8000
Па
Па.
Находим коэффициенты линейной модели
,
(3.6);
Па. (3.7)
Рассчитываем давление на выходе из
элемента, соответствующее
,
используя полученное уравнение модели
Па.
Проверяем полученное значение
,
подставляя его в поле {List
of parameter
values} при установке
параметров решателя, а затем повторяя
расчет. В результате получаем
Па,
Па.
Сравнение вычисленного значения
перепада давления с
показывает, что оно отличается на -0,09%,
что находится в пределах требуемой
точности (0,5%).
Интерпретируем полученное решение
Определяем силы, действующие на стенку канала, которая в геометрии модели соответствует границе №4. Для этого командой [Boundary integration], доступной на вкладке [Postprocessing] главной панели инструментов «FEMLAB», запускаем интегратор, и настраиваем его параметры так, как это показано на рисунке 3.22.
Рисунок 3.22 – Определение осевых составляющих удельных усилий,
действующих на стенку канала течения
В качестве подынтегрального выражения (Expression to integrate) поочередно задаем следующие переменные:
- общую удельную
силу по оси
z (
)
– “Total force per area, z component”;
- удельную силу
вязкости по
оси z (
)
– “Viscous force per area, z component”;
- общую удельную
силу по оси
r (
)
– “Total force per area, r component”;
- удельную силу
вязкости по
оси r (
)
– “Viscous force per area, r component”.
Значения интегралов появляются на информационном табло рабочего окна «FEMLAB», как показано на рисунке 3.22. Для получения значения удельной силы необходимо разделить значение интеграла на длину границы, по которой он берется, а также умножить полученный результат на -1
Н/м (3.8)
Составляющая удельной силы инерции по оси z вычисляется как разницу между соответствующими компонентами общей удельной силы и удельной силы вязкости
Н/м (3.9)
Составляющая удельной силы инерции по оси r вычисляется как разницу между соответствующими компонентами общей удельной силы, удельной силы вязкости и среднего давления
Н/м (3.10)
Среднее давление рассчитывается при помощи интегратора при задании давления – “Pressure” в качестве подынтегрального выражения (Expression to integrate), аналогично тому, как это было сделано на рисунке 3.20 для других переменных.
Результаты расчетов компонентов удельных сил сводим в таблицу 3.1.
Таблица 3.1 Силы, действующие на стенку канала течения
Из данных таблицы 3.1 следует, что основными отличиями турбулентного режима от ламинарного с точки зрения воздействия со стороны жидкости на стенку канала являются:
- изменение знака и величины касательного напряжения – при турбулентном режиме оно преимущественно действует против направления потока, но при этом имеет в 20 раз меньшую величину;
- появление составляющей силы вязкого трения по нормали к оси трубы, величина которой составляет около трети от величины касательных напряжений, действующих при ламинарном режиме течения, и в 6 раз превышает составляющую силы вязкого трения вдоль оси трубы при турбулентном режиме течения.
Таким образом, в результате сделанной работы охарактеризована подъемная сила, действующая со стороны жидкости на стенку и находящиеся на ней твердые частицы.