УДК 532.137.3
ОЦЕНИВАНИЕ МЕТОДОМ КРУТИЛЬНЫХ КОЛЕБАНИЙ ПЛОТНОСТИ НЬЮТОНОВСКОЙ СРЕДЫ ОДНОВРЕМЕННО С ВЯЗКОСТЬЮ
Г.П. Вяткин, И.В. Елюхина
Рассмотрены не обсуждаемые ранее вопросы, возникающие при одновременном оценивании вязкости и плотности ньютоновской среды в экспериментах с крутильным вискозиметром, и выполнены практические приложения предлагаемых методов параметрической идентификации.
Возможности измерения плотности р одновременно с вязкостью V в различных вискози-метрических экспериментах рассмотрены, а также детально изучены для случая внешней гидродинамической задачи, в [1]. Подобное обсуждение для крутильно-колебательного вискозиметра [2] выполнено в [3], но эти результаты не свободны от недостатков [4], затрудняющих их использование, например, для проверки согласованности вискозиметрических данных. Работоспособная методика одновременной оценки свойств к настоящему времени не реализована на практике. В [4] были изложены некоторые аспекты такого оценивания, а в настоящей работе детально рассмотрим основные из них и продемонстрируем эффективность применения разработанных методов для решения практических задач.
Математическая формулировка задачи, и в т.ч. рабочее вискозиметрическое уравнение представлены в [4]. Далее все размерные величины даны в системе СГС; под чувствительностью понимается у/у х = |(ду /дх) ■ х/у| , а под ошибками Ах параметров х - их относительные величины. Неизвестные параметры V и р определяются из условия минимума функции качества вида
I(V,р) =4%е ■ Г) + с1т ■ 1т2(Г) , (1)
где сЯе, с1т - весовые коэффициенты; Яе, 1т - действительная и мнимая части от функции Г [4].
Введем также безразмерные комплексы
А = 0,5МЯ2/ К, ^о = Я / й , ^ = 2Н / Я, (2)
где й = ч]уи0 / 2п - толщина пограничного слоя, Н - полувысота столба жидкости, К - момент инерции пустой подвесной системы, М - масса жидкости; Я - внутренний радиус цилиндра, т0 - период собственных колебаний в отсутствии среды.
Особенности одновременного оценивания вязкости и плотности
Причинами, затрудняющими решение задачи, в приведенной постановке являются овражистый характер функции (1) на плоскости (V,р) и высокая чувствительность рассчитываемых через 1т(Г) свойств к изменению периода. Тангенс угла между касательной к оси криволинейного оврага и осью плотности близок по значению к р . Заметим, что
¥р, х ~ ¥v, х ¥р, V . (3)
В вискозиметрическое уравнение плотность входит только в слагаемые, отражающие трение на торцах вискозиметра, не содержится в слагаемом, выражающем трение на боковой поверхности, и поэтому значение у/р v зависит, прежде всего, от отношения высоты к радиусу. Так, например, для длинного цилиндра овраг параллелен оси р и значение р не влияет на оценку V ; если число смачиваемых жидкостью торцов а = 1, то величина у/р v приблизительно в 2 раза
выше по сравнению с таковой для а = 2 . На рис. 1а показан поворот оси в зависимости от отношения вкладов в момент трения от торцов и боковой поверхности при базовом расчетном варианте задачи (*): Я = 1, К = 50, т0 = 5, М1 = 3, V = 0,01, р = 1, а = 2 (А = 0,03, = 11,21, X = 0,956)- кривая 1; кривая 2- а = 1, М1 = 3 ; 3 - а = 2, М2 = 30 ( А = 0,3, х = 9,56).
v-10
1.65
0.95 -
v-10
1.15 -
0.85 -
Рис. 1. Зависимость положения оси оврага от условий эксперимента
Из анализа чувствительности, проведенного в предположении независимости Ар и Av, установлено, что при cRe = 1, cIm = 0 функции yv x ~ 1.. .5 (x - K, M, R). Функции чувствительности вязкости к параметрам установки и колебаний, рассчитанные через действительную и мнимую части отличаются несущественно, кроме чувствительности к периоду, значения которой могут быть в 102 ...103 раз выше для Im(F) по сравнению с Re(F). Характер чувствительности вязкости к параметрам колебаний: периоду т и декременту S , в зависимости от условий опыта широко обсуждался в литературе с момента разработки метода (см., например, [2, 5]). Вспомним, что величина т/т0 уменьшается с ростом в основном на интервале е (2.12), а поведение S зависит от : при некотором реализуется пик (для длинного цилиндра это значение ~ 4,3, а с уменьшением высоты оно смещается в сторону больших ), т.е. в районе максимума чувствительность вязкости к декременту очень высока, также как и к ошибкам в периоде при , близких к 0 и высоких ; с уменьшением A эти чувствительности также растут.
Учитывая (3) и то, что yv т при определенных условиях эксперимента может достигать
102 ...103, получаем (при почти горизонтальном положении оврага на рис. 1а) ошибку в оценке
плотности Ар ~ 102.. ,103 % при Ат ~ 0,1 %. В случае одновременной оценки ситуация еще более
усугубляется; к тому же, помимо ошибок в параметрах колебаний в экспериментальных данных присутствуют ошибки и в параметрах установки, некоторые из которых могут составлять ~ 1%
2
(в K или в R при использовании керамического тигля), что может приводить к Ар ~ 10 %.
На рис. 1б овраги 1а и 2а соответствуют 0,99R, 1б и 2б - 1,01R , остальные параметры отвечают (*). Точки пересечения оврагов, соответствующих одинаковым радиусам (точки 1 на рис. 1б), имеют на плоскости (р, v -102) координаты: (0,98; 0,97), (1,1), (1,04; 1,02), начиная с нижней (т.е. отклонение от истинного значения не более 4%). Для AR ~ 1% при одновременном определении свойств: Ар ~ 32 %, Av ~ 15 % при M = Mj (при M2: Ар ~ 52 %, Av ~ 8 %) - этим значениям отвечают минимумы f (1) (точки 2); при независимом оценивании: Ар ~7%, Av ~ 5 % - при M = M1 (при M2 : Ар ~ 30 %, Av ~ 2,5 %) - эти точки близки к точкам пересечения осей с прямыми v = const, р = const (точки 3). При некоторой совокупности условий эксперимента и значений ошибок в измеряемых параметрах, овраги сдвигаются сильнее, и точка их пересечения может лежать дальше от истинного значения, чем минимум на оси.
Вяткин Г. П., Елюхина И. В.
Оценивание методом крутильных колебаний плотности ньютоновской среды одновременно с вязкостью
Идентификация свойств и оптимальное планирование эксперимента
Варианты оценивания плотности одновременно с вязкостью. При определении из вискози-метрической функции только вязкости среды требуется одно уравнение (Яе(Г) = 0 или 1т(Г) = 0), наиболее подходящее для этой цели в конкретной задаче. При оценке дополнительно и плотности уравнение необходимо дополнить вторым соотношением (например, традиционно, в т.ч. и в [3] проводится совместное решение для действительной и мнимой частей). Далее для оценки V и р используем следующие три варианта.
1. Предлагаемый в настоящей работе метод, когда идентификация выполняется с использованием только Яе(Г), а при планировании оптимального эксперимента в терминах теории чувствительности на множестве всех допустимых выборок измерений определяются точки, минимальное число которых (а именно, две) обеспечивает однозначное, устойчивое и надежное решение задачи. Точки отвечают различным массам образца и находятся из анализа в терминах матрицы Якоби
д5(1) д5(1) дт(1) дт(1)
др дv др дv
+ с ^ тах
д5(2) д5(2) дт(2) дт(2)
др дv др дv
(4)
где (1), (2) - номера измерений; е5 и ст - весовые коэффициенты, зависящие от точности измерения т и 5 и чувствительности V и р к ним. При выборе условий согласно (4) учитывается, что
- оси оврагов должны по возможности максимально расходиться (дно оврага становится менее пологим);
- чувствительность свойств к ошибкам в параметрах колебаний растет при уменьшении А (т.е. при уменьшении высоты и прочих равных условиях)
- величина р уменьшается с ростом х (влияние торцов ослабевает).
Элементы матриц в (4) также могут включать дополнительные весовые коэффициенты в зависимости от того, какой из параметров: V или р, должен быть оценен наиболее точно, или во избежание преобладания какой-либо из точек (так как при изменении М меняется наклон откосов оврага). В начальном приближении без анализа в рамках (4) можно принять хтах ~Юхт1п. Обратим внимание, что точка 1 с минимальной высотой, помимо перечисленного, определяется влиянием вторичных течений, возникающих у торцов, на закон колебаний. Ниже обсудим этот вопрос отдельно.
2. Метод, в котором используется только одна экспериментальная точка и оценивание производится по Яе(Г) (частный случай варианта 1 при хтах = Хть).
3. Традиционный метод оценивания (по одной точке с использованием и Яе(Г), и 1т(Г)).
Заметим, что в варианте 1 чувствительность вдоль оси оврага определяющего значения не
имеет, так как в большей степени ищутся точки пересечения оврагов, но они могут дальше отстоять от истинного значения, чем в варианте 2, основной недостаток которого - меньший наклон /(V, р) вдоль оси. При расчете по варианту 1 при Дт ~ 10-5 и Д5 ~ 10-3 (с учетом того, что точность измерения 5 обычно на 1-2 порядка ниже, чем т) и , лежащем не в окрестности
, ошибка Др < 5% (при Д5 ~10-4 - Др < 1%), что позволяет проверить согласованность вис-козиметрических данных. В варианте 2 ошибка Др при различных условиях эксперимента лежит
в широком диапазоне 1.. .10 %.
Влияние вторичных течений на колебания вискозиметра. Вторичные течения достаточно активно исследуются для внешней гидродинамической задачи (см., например, [6] и пр.), а для внутренней этот вопрос поднимался в [2, 3, 5], но представленные замечания носят только оценочный характер. Математическую модель при учете осевой и радиальной компонент скорости жидкости представим в естественных переменных скорость - давление в следующем виде:
1) уравнение движения цилиндра
с
5
d 2 а dT'
- + а = -
£оП0 015£ " «
£0
л 4 А £0П0
1 [г '
- )№^
£2 d£
П0
2) уравнения движения жидкости
дФ 5Ф ЯФ 1 д ( 5Ф
-+ V-+ W-=--1 £-
дT д£ дщ £ д£ { д£
+
д 2Ф дп2
+ 8.
(5)
(6)
где коэффициенты обобщенного уравнения (6)
>= [V, V, W]т, 8 = [-(сР/д£)- V/£2 + U2/£, - V/£2 - Ш/£, -(сР/дщ)- G]T
3) уравнение неразрывности
V дV дW п
— + — +-= 0 ;
£ д£ дщ
(7)
4) начально-краевые условия для (5)-(7) Т = 0: а
6°, = 0, V = V = W = 0; щ = 0, щ = щ0: V = W = 0, V = £ ; £ = 0: V = V = = 0; £ = £0 : V = W = 0, V = ^£0 ■
(8)
Здесь
V V ^ =
,5(р
Т = , Р =
Р
G =
£ = - п = — п = - • d d d
^0 dq0
g - ускорение свободного падения; ? - время; р - давление; q0 - циклическая частота колебаний пустого цилиндра; - - радиальная координата ( - = 0 на оси цилиндра); 5 ,5 - радиальная, азимутальная и осевая компоненты скорости; г - осевая координата ( г = 0 на нижнем торце цилиндра, г = 2Н - на верхнем); а - угловое смещение цилиндра из положения равновесия; выражение в правой части (5) - момент сил, приложенных к цилиндру со стороны среды. В (5)-(8) предполагалось, что скольжение между жидкостью и внутренней поверхностью цилиндра отсутствует, течение осесимметричное, затухание колебаний в отсутствии среды пренебрежимо мало, ньютоновская несжимаемая жидкость смачивает оба торца цилиндра. При численной реализации дискретизация уравнений проводилась методом контрольных объемов, согласование поля давления с полем скорости - с помощью 81ЫРЬБ-алгоритма, конвективные члены представлялись по рШС-схеме.
Установлено, что при малых колебаниях, реализуемых в крутильно-колебательном методе, после прохождения переходных процессов осевая и радиальная составляющие скорости не оказывают существенного влияния на закон колебаний в диапазоне величин х обычно проводимых экспериментов (см., например, [7, 8] и пр.). При высокой точности измерения периода и декремента затухания (не хуже 10-5 ...10-3) целесообразно ограничиться хшт >~ 0,7...0,9, а при необходимости работы с вискозиметрами малой высоты с х <~ 0,1...0,2 следует при заданных условиях эксперимента найти оценки коэффициентов для коррекции параметров колебаний с учетом течений, возникающих у торцов. Полученные результаты соответствуют распространенным для режима затухающих колебаний значениям £0 ~ 8,5...25 ( А = 0,02...0,2).
Для сильновязкого приближения [2] при прочих равных условиях вторичные течения развиты в большем объеме вискозиметра, а с ростом £0 значительно усиливается интенсивность таких течений вблизи тигля, что искажает поле азимутальной скорости у поверхностей, изменяя момент сил трения в (5) и, следовательно, зависимость а = а(Т) по сравнению с аналитическим
решением. Кроме того, как уже указывалось выше, в определенном диапазоне значений £ 0 наблюдается высокая чувствительность свойств к ошибкам в параметрах. В этих случаях граница значений хш1п, обеспечивающая работоспособность методики оценивания свойств жидкости, смещается вверх. При планировании оптимального эксперимента по выбору минимально допустимой высоты тигля хшп учитывались также длительность и интенсивность выхода на регулярный режим колебаний, определяемые при заданных £ 0 и щ0 комплексом А. Так, при больших
Вяткин Г.П., Елюхина И.В. Оценивание методом крутильных колебаний плотности
ньютоновской среды одновременно с вязкостью
5 переходные процессы могут занимать значительную часть временного интервала, на котором возможна корректная регистрация а , а при низких 5 после прохождения таких процессов реализуются высокие амплитудные значения углового смещения (а ~0,1) и отнесенные к азимутальной компоненте значения осевой и радиальной составляющих скорости, что сильнее проявляется с ростом % 0.
Практические приложения разработанных методов
Пример 1. Примем размеры тигля, согласно [3] не удовлетворительные для одновременной оценки V и р, а точность измерения т и 5 хуже, чем в [3]; к тому же, пусть значение А будет
мало (что затрудняет даже оценку V), а % 0 близко к % 05 . Положим: Я = 0,5, МтП = 0,4, К = 15, т0 = 5, V = 0,01, р0 = 1 (А = 0,0033, %0 = 5,6, хтт = 1,019). Проведем учет ошибок только т и 5 , как и в [3] (Дт = Д5 = 0,001).
На рис. 2 приведены зависимости ошибок Др = Др (М) (здесь и далее
Др = |(р - р0)/ р0| -100%); номера кривых соответствуют номерам вариантов, в варианте 1 -М = Мтах ; расчет выполнен для наиболее распространенного случая а = 1. Значения Др на
рис. 2 - максимальные из получаемых при всех возможных комбинациях т ± Дт , 5 ± Д5 . При расчете по традиционному варианту 3 при определенном соотношении знаков в ошибках Дт, Д5 оказывается р ~ 0 ; значение V ~ 10-4 (при М = МтП ), с ростом высоты влияние р на величину вязкости падает, и при М = 4 величина Дv ~ 20 %. Видно, что наилучшие значения достигаются в варианте 2 при малой высоте, но при учете ошибок во всех входных данных и высокой точности измерения декремента этот метод обычно менее эффективен предлагаемого Рис-2 Зависимость °шиб°к в пл°™°сти
, от высоты заполнения тигля
варианта 1.
Пример 2. Выполним обработку реальных экспериментальных данных для воды ( р0 = 1). Здесь ДЯ ~ 0,005, ДК ~0,02, ДМ ~ 0,001, Дт ~ 0,0005, Д5 ~0,02; радиус тигля Я = 1,525, высота хтах = 3,2. Обозначим номера расчетных точек как, например, 23 - соответствующая второму сверху значению К и массе М3 (см. таблицу). При одновременном определении V и р традиционным методом в расчетах 22, 32, 41, 42, 51, 52 значения р ~ 0 , в 11, 13, 21, 43 - р > 2 , а в 12 и 53 - Др < 10%. В случае варианта 2 значения р ~ 0,88 ± 0,1. При расчете же по варианту 1 построчно с использованием крайних значений масс М1 и М3 в среднем Др ~ 20%, а для второй строки даже при выборе оптимальных высот в рамках (4) р ~ 0 . Итак, предлагаемый метод оказался хуже, чем традиционный. Проанализируем, случайность ли это.
Обратимся к теории чувствительности и определим, какие максимальные ошибки могут возникать в вариантах 1-3. При использовании только действительных частей в (1) для указанного диапазона х и К значения Др ~ 20 ± 5%, а при Д5 < 0,01, ДК < 0,01, ДЯ < 0,001 может оказаться Др < 5 %. Заметим, что полученные оценки р укладываются в этот интервал. В варианте 1 добавляются еще дополнительно ошибки в т, 5, М для второй точки. Ошибка здесь повышается также и из-за роста у/р 5 (ошибка от Д5 в плотности Др ~ 30 %; в варианте 2 - Др ~ 7 %), так как определение свойств среды по ) учитывает, прежде всего, значения 5 и Д5 . В этом методе желательно точнее измерять декремент: при Д5 < 0,001 и прочих равных условиях можно надежно получить Др < 5 %.
Таблица
K M = 0 M1 = 15,115 (X = 1,357) M2 = 20,197 (X = 1,813) M3 = 25,259 (X = 2,267)
Т0 ^0 Т 5 A Т 5 A Т 5 A
118 7,131 0,0216 14,3 7,263 0,101 0,149 7,292 0,133 0,199 7,316 0,149 0,249
146 7,922 0,0208 13,6 8,045 0,086 0,12 8,053 0,111 0,161 8,088 0,13 0,201
196 9,158 0,0138 12,6 9,294 0,072 0,089 9,233 0,0904 0,12 9,271 0,106 0,15
236 9,998 0,0126 12,1 10,078 0,0639 0,074 10,09 0,0792 0,099 10,22 0,097 0,124
274 10,773 0,0137 11,6 10,853 0,0578 0,064 10,85 0,0726 0,086 10,92 0,084 0,107
В примере наблюдается ситуация иная, чем на рис. 2: точки пересечения оврага расположены дальше от истинного значения, чем минимум Re(F) на оси оврага. Традиционный вариант здесь совершенно не применим к одновременной оценке свойств: помимо достаточно высоких значений Ар от ошибок в K и R (~ 15 %) ошибки Дг ведут или к р ~ 0 , или к Др > 500%; при
оценке только вязкости Av ~ 5%. Итак, ошибки Аx оказали такое влияние (т.е. некоторые из них увеличили значения вязкости и плотности, а некоторые уменьшили), что даже традиционный вариант оказался лучше предлагаемого (например, точка 12 и пр.); эти ошибки легко прогнозируемы и контролируемы при математическом моделировании условий крутильных колебаний.
Заключение
Итак, в настоящей работе показано, что
1) в зависимости от ошибок наблюдаемых параметров и условий эксперимента лучшим может быть вариант идентификации по действительной части вискозиметрического уравнения с использованием одной или двух точек и для проверки согласованности следует выбрать именно его, а затем путем моделирования оценить доверительные интервалы для оценок свойств в этих случаях и посмотреть, укладываются ли в них найденные из эксперимента значения v и р ;
2) для исключения эффектов вторичных течений и возможности использования традиционных расчетных выражений эксперименты по проверке согласованности вискозиметрических данных рекомендуется проводить в области слабовязкого приближения при высотах тигля, не меньших радиуса, и в предлагаемом методе идентификации свойств принять Xmn ~ 1.
Авторы признательны Бескачко В.П. за предоставление обсужденных выше опытных данных.
Литература
1. Krall A.H., Sengers J.V. Simultaneous measurement of viscosity and density with an oscillating-disk instrument// Int. J. of Thermophys. - 2003. - V. 24. - № 2. - Р. 337-359.
2. Швидковский Е.Г. Некоторые вопросы вязкости расплавленных металлов. - М.: ГИТТЛ, 1955.-206 с.
3. Nieuwoudt J.C., Sengers J.V., Kestin J. On the theory of oscillating-cup viscometers// Physica А. - 1988.-V. 149.- Р. 107-122.
4. Елюхина И.В. Планирование оптимального эксперимента по одновременному определению вязкости и плотности ньютоновской среды// Вестник ЮУрГУ. Серия «Математика, физика, химия». - 2001. - Вып. 1. - № 7(07) - С. 71-76.
5. Grouvel J.M., Kestin J. Working equations for the oscillating-cup viscometer// Appl. Sci. Res. -1978.-V. 34.- Р. 427-443.
6. Folse R.F. Observations of secondary flows generated by a torsionally oscillating sphere in linearly stratified fluids// Phys. Fluids. - 1994. - V. 6. - Р. 537-540.
7. Vollmann J., Riedel D. The viscosity of liquid Bi-Ga alloys// J. Phys.: Condens. Matter. - 1996. -V. 8.- Р. 6175-6184.
8. Wang D., Overfelt R.A. Oscillating rap viscosity measurements of aluminum alloys: A201, A319 and A356// Int. J. of Thermophys. - 2002. - V. 23. - № 4. - P. 1063-1076.
Поступила в редакцию 27 октября 2005 г.