Друзья, поддержите, если возможно советом. За последние неск. дней все мозги сломал, но так и не нашёл ошибки. Подозреваю, что банально замылился глаз и не вижу того, что лежит под носом.
Итак, есть стандарт. В нём приведены чёткие уравнения и алгоритмы расчёта плотности нефтепродуктов, в зависимости от температуры/давления. Вот примерная схема расчётов:
Скрытый текст
С этим вроде как всё ясно и понятно. Изначальная задача: вывести р15 из замеров топлива, находящегося в произвольном состоянии. На основании вышепривёдённых данных, родился следующий текст в VBA:
И вроде бы всё хорошо. Но, не всё хорошо. Считает не правильно. В том же стандарте приводится пример для расчёта. Следуя которому - мною написанное даёт иные результаты. Вот пример:
Скрытый текст
Прицепил файлик, в котором играюсь. Буду крайне признателен, если ткнёте носом в ошибку!
panix1111 написал: округляете же сами, типы переменных для такой точности надо указывать
Проблема не в округлениях. Там точность в худшем случае пострадает. Проблема в том, что он в принципе неверно считает.
Грубо говоря, при пересчёте из более высокой температуры в более низкую - плотность должна возрастать (что, как минимум, логично), а она - уменьшается. Другими словами: если при температуре +27, плотность солярки 840, то при +15 она должна быть +845, а не - 835.
Если взять конкретный пример и первого поста: там при t=27,3 плотность = 836,15 в результате первого приближения, при приведении к t=15 - формула должна выдать - 843,62, а она выдаёт - 827.31.... в чём собственно и состоит затык...
Изменено: Святослав И - 29.10.2017 13:47:52(уточнение ответа)
Пошаговый прогон показывает, что проблема в ф-ле Yt, так b15 выдает значение как и в ГОСТ-примере, а вот Yt = 8.1328 (должно быть 8,148) Хотя, вроде б ф-ла записана правильно. Но ответ то не сходится, а потом и всё остальное будет неверно. Понял, чего точно не хватет - второй температуры, первая при которой меряли д ля расчета приближения, а вторая(для которой расчет и идет) должна использоваться на последующих шагах
У меня возникло такое ощущение, что в стандарте где-то что-то пропущено. У меня в результате вычислений получаются те же цифры, что и у Вас. При этом, более-менее близкие к примеру цифры получаются только при убирании знака минус в формуле 1 у коэффициента Б15.
Программный код, как и яды, лучше тестировать по капельке Люблю изобретать велосипеды с колесами произвольной формы
panix1111 написал: Пошаговый прогон показывает, что проблема в ф-ле Yt, так b15 выдает значение как и в ГОСТ-примере, а вот Yt = 8.1328 (должно быть 8,148)
В прицепленом файлике, на листе "Руками" - я простыми формулами посчитал - и там, в первом приближении - значение коэффициента Yt - в точности совпадает со значением из примера.
Цитата
panix1111 написал: Понял, чего точно не хватет - второй температуры, первая при которой меряли д ля расчета приближения, а вторая(для которой расчет и идет) должна использоваться на последующих шагах
Хм... Так вроде в алгоритме указано использовать именно замеренную температуру то есть PtP... Пошёл читать в сто первый раз алгоритм...
Вот же:
Цитата
Примечания При определении значения плотности p15, в формулы (1) и (3) подставляют значения температуры и давления нефти или нефтепродукта, при которых была измерена плотность PtP НО При расчете значения плотности PtP в формулы (1) и (3) подставляют значения температуры и давления нефти, при которых требуется определить плотность.
В рассматриваемом примере мы выводим именно p15
Цитата
Hypohelix написал: У меня возникло такое ощущение, что в стандарте где-то что-то пропущено. У меня в результате вычислений получаются те же цифры, что и у Вас. При этом, более-менее близкие к примеру цифры получаются только при убирании знака минус в формуле 1 у коэффициента Б15.
Угу.... тоже игрался с этими минусами... И потом - вряд ли в госте напортачили... Я сверялся с предыдущими гостами - там практически те же самые формулы используются...
Function P_to_15(dens As Double, temp_start As Single, Optional p As Single = 0) As Double
Dim old_dens As Double, b15 As Double, Yt As Double, pTp As Double
Const k0 = 613.9723
Const k1 = 0
Const k2 = 0
old_dens = 0
pTp = dens
t = temp_start
Do Until pTp - old_dens <= 0.01
b15 = ((k0 + k1 * pTp) / pTp ^ 2) + k2
Yt = 10 ^ -3 * Exp(-1.6208 + (0.00021592 * t) + _
((0.87096 * 10 ^ 6) / pTp ^ 2) + ((4.2092 * t * 10 ^ 3) / pTp ^ 2))
old_dens = pTp
'Exp(b15 * (t - 15) * Abs((1 + 0.8 * b15 * (t - 15)))))
pTp = (old_dens * Exp(-b15 * (t - 15) * Abs((1 + 0.8 * b15 * (t - 15))))) / (1 - Yt * p)
Loop
P_to_15 = pTp
End Function
Да, разобрался, таки надо указывать давление, если ставить ссылочку на третий параметр, то выходит так же как и в вашем варианте "вручную" = 828.7441072774
panix1111 написал: что то там не то с вычислениями
Угу, мягко говоря - не то... Настолько не то, что оно уходит в противоположном направлении. С другой стороны, вот ещё раз формула:
как вот её ещё можно понять иначе? У меня, грешным делом, была мысль, что скобки квадратные-фигурные каким-то образом могут влиять на очерёдность вычислений...
Святослав И написал: скобки квадратные-фигурные каким-то образом могут влиять на очерёдность вычислений
квадратные - это модуль, я его в последнем варианте на всяк случай добавил в ф-лу (Abs), но всё равно что-то не сходится с расчетами из примера. Что могу подсказать - есть готовые таблицы для расчета, может их как-то приспособить?
panix1111 написал: есть готовые таблицы для расчета, может их как-то приспособить?
Да, такой вариант, за неимением лучшего, я уже приспособил и что-то оно вычисляет. Но. Погрешности к сож., слишком велики. Велики на столько, что всю ценность вычислений сводят на "нет". Ну то есть, грубо говоря, на бензовозе в котором плещется 22-23 тонны солярки - литров 200 свободно "гуляет" туда/сюда только на этих знаках после запятой... Соответственно, какие-то недоливы/сливы можно даже не пытаться ловить, соответственно - про какой-либо вменяемый учёт/контроль можно забывать... Просто я покурил форумы метрологов и там вполне чётко прослеживается мысль, что только вот эти формулы...
При определении значения плотности p15 в формулы (1) и (3) подставляют значения температуры и давления нефти или нефтепродукта, при которых была измерена плотностьPtP. При расчете значения плотности PtP в формулы (1) и (3) подставляют значения температуры и давления нефти, при которых требуется определить плотность.
panix1111, так я ровно о том же самом написал в #7-ом сообщении. Фишка в том, что мы как раз вычисляем именно р15, а не РtP... В алгоритме подстановок вроде об этом явно написано, или я как-то не так его читаю?
Цитата
3) Значение P15, вычисленное в первом приближении, подставляют в формулы (2) и (3) и вычисляют во втором приближении значения B15 и Yt, соответственно
Хотя вот начиная с этого пункта возникает некоторая неопределённость. Но опять же. Начиная с самой первой подстановки полученных b15 и Yt (значения которых, кстати, вычисляются чётко по примеру) - формула возвращает неверное значение, т.е. уже на самом первом шаге.
Function Beta15(D15 As Double, K0 As Double, K1 As Double, K2 As Double)
Rem Возвращает Коэффициент объемного расширения при t=15
Rem D15 - плотность при температуре 15 и избыточном давлении = 0
Beta15 = (K0 / D15 + K1) / D15 + K2
End Function
Function GammaT(D15 As Double, T As Double)
Rem Возвращает Коэффициент сжимаемости при температуре T
Rem D15 - плотность при температуре 15 и избыточном давлении = 0
GammaT = 10 ^ -3 * Exp(-1.6208 + 0.00021592 * T + (0.87096 * 10 ^ 6 + 4.2092 * T * 10 ^ 3) / D15 ^ 2)
End Function
Function Density15(DTP As Double, T As Double, P As Double, K0 As Double, K1 As Double, K2 As Double)
Rem DTP - известная, измеренная прибором (ареометром) при температуре T и давлении P, плотность нефтепродукта
Rem Рассчитывает плотность при температуре 15 и избыточном давлении = 0
Dim D15 As Double, B15 As Double, DTemp As Double, ICnt As Long
D15 = DTP
Do
DTemp = D15
B15 = Beta15(D15, K0, K1, K2)
D15 = DTP * (1 - GammaT(D15, T) * P) * Exp(B15 * (T - 15) * (1 + 0.8 * B15 * (T - 15)))
ICnt = ICnt + 1
If ICnt > 100 Then
Density15 = CVErr(xlErrNA): Exit Function
End If
Loop Until Abs(D15 - DTemp) < 0.01
Density15 = D15
End Function
Function DensityTP(D15 As Double, T As Double, P As Double, K0 As Double, K1 As Double, K2 As Double)
Rem Рассчитывает плотность нефтепродукта при температуре T и избыточном давлении P
Rem D15 - плотность при температуре 15 и избыточном давлении = 0
Dim B15 As Double
B15 = Beta15(D15, K0, K1, K2)
DensityTP = D15 * Exp(-B15 * (T - 15) * (1 + 0.8 * B15 * (T - 15))) / (1 - GammaT(D15, T) * P)
End Function
Нас нае... обманули. Расходимся. Оказывается (!!!) для того чтобы формула из стандарта, вот эта:
Работала правильно для приведения PtP в P15, её всего лишь (!!!) нужно перевернуть. Вот так:
И тогда всё внезапно начинает правильно работать. То, что в ГОСТе забыли об этом написать, ну так то мелочи. Нет слов.
Всем спасибо за посильный вклад!
P.S. Это мне народ на метрологическом форуме подсказал. Для них это - нечто само собой разумеющееся. P.P.S. ...Как и здесь для отдельных людей - другие вещи.
Святослав И, извините, а что Вы рассчитывали на листе Анализ ? Плотность при t=15 и P=0 или какую ?
Ещё, в РЕКОМЕНДАЦИИ написали: "Значение плотности ρ15 находят методом последовательных приближений, используя итерационный метод прямых подстановок", но забыли упомянуть, что ρ15, в этом случае определяется по "обратной" формуле:
Т.е. ниже Вы пытаетесь определить плотность опять при тех же 27.3°С, 2.45 МПа, и измеренной плотности 836,15 кг/м3 ???
В файле, то что вы VBA проделываете с помощью цикла, я просто разложил на листе - для большей наглядности. Так лучше видно - на какой итерации какие результаты выходят. Собственно, debug.print в консоль (в вашей функции) приносит ровно те же результаты. То, что я пытался делать - привести плотность при 27,3°С в плотность при 15°С, при том же P (т.е. 2,45 МПа) для пересчёта в другое значение P - там есть свои формулы. При этом чётко следовал, рекомендациям, о том, что дескать нужные значения просто подставить в ту же формулу (1). То, что вычислять нужно на перевёрнутой формуле - ребята забыли упомянуть (равно как и привести её). Ну что ж, бывает...
Ещё раз, всем спасибо за посильное участие и помощь.