Раздел 1. Основные типы лечебно-профилактических учреждений 4 страница

После выбора типа решателя опцией {Solver}, задаем его основные свойства на вкладке [General], как это показано на рисунке 2.4. Для типа решателя “Time Dependent” на этой вкладке устанавливается опция {Time stepping} (шаги по времени), для которой параметром {Times:} задаются временные слои в формате начальный момент:шаг:конечный момент, а также параметрами {Relative tolerance:} и {Absolute tolerance:} определяется относительная и абсолютная ошибки на каждом временном шаге. Установка галочки в поле {Allow complex numbers} позволяет применять в расчетах комплексные числа, что необходимо при задании коэффициентов дифференциальных уравнений в комплексном виде.

На вкладке [Time Stepping] параметром {Times to store in output:}устанавливаются моменты времени, для которых будет выводиться решение, в том числе настраивается их автоматический выбор самим решателем. Значение этого параметра по умолчанию – “Specified Times”, то есть будут выведены решения для моментов времени, определённых на вкладке General. Для того чтобы увидеть результаты вычислений для всех шагов решающего устройства надо выбрать “Time steps from solver”. Вообще, решающее устройство выбирает шаги самостоятельно в зависимости от динамики решения, то есть игнорирует их значения, заданные на вкладке [General]. Чтобы решающее устройство учитывало этот список, например, если внешние воздействия импульсные и решатель может «проскочить» мимо них, необходимо установить значение параметра {Time steps taken by solver:} “Strict” – тогда решатель будет использовать только эти шаги, или “Intermediate” – решатель будет использовать как свободные шаги, так и шаги, перечисленные на вкладке General. Значение этого параметра по умолчанию “Free”, т.е. допускаются только свободные шаги. Если надо принудительно задать шаг по времени нужного размера, то это делается посредством опции {Manual tuning of step size}.

Для типа решателя “Eigenvalue” необходимо задать количество собственных чисел уравнения параметром {Desired number of eigenvalues:}, а также параметром {Search for eigenvalues around:}определяется область их поиска.

Для нелинейных решателей появляется вкладка [Nonlinear], где можно указать количество итераций. Для сильно нелинейных задач рекомендуется поставить галочку напротив параметра {Highly nonlinear problem} и увеличить количество итераций сверх значения 25, принятого по умолчанию.

Для всех режимов, кроме режима «Time Dependent», можно установить галочку напротив параметра {Adaption}, что позволяет уточнить сетку на основе уже полученного решения. Это рекомендуется делать, если физика и геометрия достаточно сложна, и не очень ясно как задавать параметры сетки. Однако время вычислений с адаптивной сеткой больше, нежели с простой.

Если расчетная матрица симметрична, то можно ускорить расчет, установив опцию {Matrix symmetry} в значение “Symmetric”. Чаще всего матрицы симметричны в задчах чистой кондукции или диффузии. Матрицы заведомо несимметричны в мультифизических моделях, моделях связанных с конвекцией или уравнением Навье-Стокса.

Какая бы задача не рассматривалась, при решении она сводится к системам линейных уравнений. Большая часть времени расчета уходит на решение этих систем, чем занимается вычислитель корней систем линейных уравнений, являющийся частью решателя «FEMLAB». Конкретный вид вычислителя определяется опцией {Linear system solver}, которая может принимать следующие значения: а) “Direct (UMFPACK)”; б) “Direct (SPOOLES)”; в) “Cholesky (TAUCS)”; г) “Conjugate Gradients”. По умолчанию стоит “Direct (UMFPACK)”, который является наиболее устойчивым, но отнимает много ресурсов компьютера, поэтому для моделей требующих длительного расчета стоит попытаться воспользоваться альтернативным решателем. В частности, если решатель “Direct (UMFPACK)” не работает или работает недопустимо долго, имеет смысл задействовать решатель “Direct (SPOOLES)”, так как он требует меньше оперативной памяти, однако при этом могут появиться проблемы со стабильностью решения. В крайнем случае, выбирается итеративный решатель “GMRES”, характеризующийся наихудшей сходимостью. Для положительно определенных систем с симметричными матрицами оптимальный результат дает решатель “Direct Cholesky (TAUCS)” или итеративный “Conjugate Gradients”.

3 ПРАКТИКУМ

3.1 Расчет ламинарного установившегося течения в трубах

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

Рассмотрим этапы расчета, применительно к промывке биополимерным соленасыщенным буровым раствором (БПС-БР) плотностью r_пж =1180 кг/м³, который прокачивается по бурильным трубам марки ПН-127´9 общей длиной 1000 м с расходом Q =0,023 м³/с.

Путем статистической обработки результатов реометрии БПС-БР получаем реологическую модель Carreau в следующем виде

. (3.1)

Работа с «FEMLAB» начинается с установки глобальных параметров модели. К ним относятся:

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

- наименование модуля, подраздела модуля и модели;

- тип расчета;

- вид конечных элементов.

Глобальные параметры модели устанавливаются в окне “ Model Navigator ”, которое открывается автоматически при запуске «FEMLAB». Первоначально в этом окне необходимо установить размерность задачи, так как от нее зависит перечень доступных для расчета моделей. Кроме того, в навигаторе моделей задается размерность, в данном случае "Axial symmetry (2D)", выбирается стационарный поток неньютоновской жидкости в модуле " Chemical engineering module ", и указывается вид конечного элемента "Langrange p3-p2", как показано на рисунке 3.1.

Рисунок 3.1 – Окно навигатора моделей

Задаем геометрию элемента (рисунок 3.2) в режиме "Draw mode", в который «FEMLAB» переходит автоматически после выбора модели. В соответствии с заданной размерностью, геометрией является половина продольного сечения внутреннего объема элемента трубы, представляющая собой прямоугольник.

Рисунок 3.2 – Окно редактора геометрии (Draw mode)

Для построения прямоугольника используем инструмент [Rectangle], расположенный на панели инструментов, вводя параметры фигуры в окне ее свойств, как это показано на рисунке 3.2.

Задаем свойства сплошной среды, как это показано на рисунке 3.3, перейдя в режим "Subdomain Mode", для чего выбираем команду [Subdomain settings...] на вкладке [Physics] панели инcтрументов «FEMLAB».

Рисунок 3.3 – Задание свойств жидкости в окне "Subdomain settings"

Задаем граничные условия, как это показано на рисунке 3.4, перейдя в режим "Boundary Mode", для чего выбираем команду [Boundary settings...] на вкладке [Physics] панели инструментов «FEMLAB».

Рисунок 3.4 – Задание граничных условий в окне "Boundary settings"

Строим сетку, используя кнопки <Initialize mesh> и <Refine mesh> на панели инcтрументов «FEMLAB». В результате получается сетка, состоящая из 2176 элементов, показанная на рисунке 3.5.

Рисунок 3.5 – Конечноэлементная сетка, отображаемая в режиме "Mesh Mode"

По самому короткому ребру геометрии насчитывается 15 элементов сетки, что соответствует умолчаниям программы и минимальным требованиям к разбиению.

Устанавливаем параметры решателя для получения начального приближения, как это показано на рисунке 3.6, используя кнопку <Solver parameters> на панели инструментов «FEMLAB».

Рисунок 3.6 – Задание параметров решателя в окне "Solver Parameters"

Значение давления на входе в элемент в поле {List of parameter values} задаем произвольно, так как – это лишь начальное приближение. Решатель выбираем линейный параметрический, так как иначе решение вероятнее всего не сойдется. Поскольку задача решается в 2-D размерности, что не требует много времени и ресурсов, выбираем наиболее устойчивый прямой решатель UMFPACK.

Открываем менеджер решения, используя кнопку <Solver manager> на панели инcтрументов «FEMLAB», и устанавливаем опции для получения начального приближения, как показано на рисунке 3.7.

Рисунок 3.7 – Установка опций менеджера решения для получения

начального приближения

После установки опций не забываем нажать кнопку <Apply> в нижней части окна "Solver Manager" для того, чтобы внесенные изменения вступили в силу.

Запускаем расчет, используя кнопку <Solve> в нижней части окна "Solver Manager". Ход расчета контролируем по данным в окне показанном на рисунке 3.8.

Рисунок 3.8 – Отображение хода решения

Процесс решения идет нормально если значение параметра “Convergence” для нелинейного решателя закономерно уменьшается по мере приближения значения параметра “Progress” к 100%. Если значение параметра “Convergence” скачком увеличивается, а параметр “Progress” при этом уменьшается вместо того, чтобы возрастать, то скорее всего модель настроена плохо и ее необходимо скорректировать, не дожидаясь достижения максимального числа итераций.

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

Рисунок 3.9 – Установка опций менеджера решения для получения

точного решения

Чтобы обеспечить возможность возврата к начальному приближению при неудачной попытке отыскания точного решения до запуска этого процесса нажимаем кнопку <Store Solution> в нижней части окна "Solver Manager", и присваиваем имя сохраняемому решению.

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

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

Рисунок 3.10 – Задание параметров решателя для отыскания точного решения

Работа с нелинейным решателем требует применения искусственной стабилизации решения. Для ее включения заходим в окно "Subdomain settings", где нажимаем кнопку <Artificial Diffusion> и в открывшемся окне устанавливаем опции, как показано на рисунке 3.11.

Рисунок 3.11 – Настройка искусственной стабилизации решения

Если решение после этого не сошлось имеет смысл попробовать на первом шаге увеличить значение параметра настройки (Tuning parameter) до 0,1-0,2, а при последующих решениях постепенно его уменьшить.

После того как точное решение получено определяем значение расхода промывочной жидкости, соответствующие перепаду давления на торцах элемента. Для этого используем команду [Boundary integration], доступную на вкладке [Postprocessing] панели инcтрументов «FEMLAB». Опции интегратора устанавливаем так, как это показано на рисунке 3.12.

Рисунок 3.12 – Задание параметров интегратора в окне "Boundary Integration"

Значение расхода 0,015972 после запуска интегратора появляется в области сообщений внизу рабочего окна «FEMLAB», где также указывается вид интеграла (Surface integral), переменная, по которой выполнено интегрирование (Expression) и пределы интегрирования (Boundary). Интегрировали в данном случае компоненту скорости жидкости, направленную вдоль оси потока, в пределах входного сечения элемента.

Поскольку полученное значение расхода меньше требуемого 0,023 м³/с требуется увеличить перепад давления на элементе. Ориентировочно для достижения нужного расхода перепад должен быть увеличен в 1,5-2 раза. Чтобы уменьшить число итераций при поиске необходимого перепада воспользуемся моделью квадратичного полинома, позволяющего найти достаточно точное приближение после получения первых четырех точек зависимости DP = f (Q). Для построения полинома выполняем расчет при следующих значениях давлениях на входе в элемент: а) 51; б) 80; в) 87,8; г) 100. Для каждого значения давления определяем расход жидкости. Полученные таким образом точки переносим в MS Excel, строим график DP = f (Q), и добавляем на нем линию тренда соответствующего типа. Результат этих построений показан на рисунке 3.13.

Рисунок 3.13 – Зависимость перепада давления в расчетном элементе от расхода

через него промывочной жидкости

Подставляя в уравнение квадратичного полинома требуемое значение расхода, вычисляем перепад давления равный 88,044 Па. Подставляем это значение в программу и убеждаемся, что расход соответствует заданному. Тогда искомые потери давления на всей длине труб составят 0,88 МПа. Результат расчета показан на рисунке 3.14.

Рисунок 3.14 – Результат решения задачи

3.2 Расчет турбулентного течения в трубах

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

Для решения задачи используются следующие исходные данные

Расчетным элементом является канал течения промывочной жидкости длиной , представляющий собой внутреннюю полость бурильной трубы типа ПН-127´10.

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

- верхняя граница средней скорости сдвига, соответствующая началу турбулентного течения – =240 с^-1;

- нижняя граница средней скорости сдвига, соответствующая началу пластического течения – =16 с^-1;

- пластическая вязкость – =0,4779 Па*с;

- реологическая модель Carreau, построенная по точкам реологического профиля, лежащим в области ламинарного течения, и пластической вязкости

;

- модель удельного повышения потерь давления при турбулентном течении

.

Плотность промывочной жидкости составляет r_пж =1120 кг/м³, ее расход при промывке скважины – Q =0,056 м³/с.

Решение

Определяем удельные потери давления, возникающие при ламинарном течении бурового раствора внутри труб, используя методику предыдущей работы. В результате вычислений получаем: а) удельные потери давления, равные =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)