Страницы: 1 2 След.
RSS
Расчёт плотности нефтепродуктов, по стандарту ГОСТ Р 50.2.076-2010
 
Здравствуйте,

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

Итак, есть стандарт. В нём приведены чёткие уравнения и алгоритмы расчёта плотности нефтепродуктов, в зависимости от температуры/давления. Вот примерная схема расчётов:
Скрытый текст

С этим вроде как всё ясно и понятно. Изначальная задача:  вывести р15 из замеров топлива, находящегося в произвольном состоянии. На основании вышепривёдённых данных, родился следующий текст в VBA:
Код
Function P_to_15(dens, temp, Optional p)
k0 = 613.9723
k1 = 0
k2 = 0
If IsMissing(p) = True Then
    p = 0
End If

old_dens = 0
new_dens = dens

Do Until Abs(new_dens - old_dens) = 0
    b15 = ((k0 + k1 * new_dens) / new_dens ^ 2) + k2
    Yt = 10 ^ -3 * Exp(-1.6208 + 0.00021592 * temp + ((0.87096 * 10 ^ 6) / new_dens ^ 2) + ((4.2092 * temp * 10 ^ 3) / new_dens ^ 2))
    old_dens = new_dens
    new_dens = (dens * Exp(-b15 * (temp - 15) * (1 + 0.8 * b15 * (temp - 15)))) / (1 - Yt * p)
Loop
P_to_15 = Application.Round(new_dens, 2)
End Function
И вроде бы всё хорошо. Но, не всё хорошо. Считает не правильно. В том же стандарте приводится пример для расчёта. Следуя которому - мною написанное даёт иные результаты. Вот пример:
Скрытый текст

Прицепил файлик, в котором играюсь. Буду крайне признателен, если ткнёте носом в ошибку!
Изменено: Святослав И - 29.10.2017 15:45:14
 
округляете же сами, типы переменных для такой точности надо указывать
Код
Function P_to_15(dens As Double, temp As Integer, Optional p As Single = 0) As Double
Dim k0 As Double, k1 As Double, k2 As Double, old_dens As Double
Dim b15 As Double, Yt As Double, new_dens As Double

k0 = 613.9723
k1 = 0
k2 = 0

old_dens = 0
new_dens = dens

Do Until new_dens - old_dens <= 0.01
    b15 = ((k0 + k1 * new_dens) / new_dens ^ 2) + k2
    Yt = 10 ^ -3 * Exp(-1.6208 + 0.00021592 * temp + ((0.87096 * 10 ^ 6) / new_dens ^ 2) + ((4.2092 * temp * 10 ^ 3) / new_dens ^ 2))
    old_dens = new_dens
    new_dens = (dens * Exp(-b15 * (temp - 15) * (1 + 0.8 * b15 * (temp - 15)))) / (1 - Yt * p)
Loop
'P_to_15 = Application.Round(new_dens, 2)
P_to_15 = new_dens
End Function

по примеру дает ответ 827.3113441474

стоп, там еще есть ошибки. плюс не учтены все коэфициенты в Юзерфкнкции
Изменено: panix1111 - 29.10.2017 13:41:33
Мы в Екселе не работаем, мы в нём живём!
 
Цитата
panix1111 написал:
округляете же сами, типы переменных для такой точности надо указывать
Проблема не в округлениях. Там точность в худшем случае пострадает. Проблема в том, что он в принципе неверно считает.

Грубо говоря, при пересчёте из более высокой температуры в более низкую - плотность должна возрастать (что, как минимум, логично), а она - уменьшается. Другими словами: если при температуре +27, плотность солярки 840, то при +15 она должна быть +845, а не - 835.

Если взять конкретный пример и первого поста: там при t=27,3 плотность = 836,15
в результате первого приближения, при приведении к t=15 - формула должна выдать - 843,62, а она выдаёт - 827.31.... в чём собственно и состоит затык...
Изменено: Святослав И - 29.10.2017 13:47:52 (уточнение ответа)
 
Цитата
panix1111 написал:
стоп, там еще есть ошибки. плюс не учтены все коэфициенты в Юзерфкнкции
Да вроде всё учёл.... всего-то три уравнения задействованы... что до ошибок - разумеется они есть ) если бы их не было и вопрос бы не стоял ))
 
Пошаговый прогон показывает, что проблема в ф-ле Yt, так b15 выдает значение как и в ГОСТ-примере, а вот Yt = 8.1328 (должно быть 8,148)
Хотя, вроде б ф-ла записана правильно. Но ответ то не сходится, а потом и всё остальное будет неверно.
Понял, чего точно не хватет - второй температуры, первая при которой меряли д ля расчета приближения, а вторая(для которой расчет и идет) должна использоваться на последующих шагах
Изменено: panix1111 - 29.10.2017 15:00:01
Мы в Екселе не работаем, мы в нём живём!
 
У меня возникло такое ощущение, что в стандарте где-то что-то пропущено. У меня в результате вычислений получаются те же цифры, что и у Вас. При этом, более-менее близкие к примеру цифры получаются только при убирании знака минус в формуле 1 у коэффициента Б15.  
Программный код, как и яды, лучше тестировать по капельке
Люблю изобретать велосипеды с колесами произвольной формы
 
Цитата
panix1111 написал: Пошаговый прогон показывает, что проблема в ф-ле Yt, так b15 выдает значение как и в ГОСТ-примере, а вот Yt = 8.1328 (должно быть 8,148)
В прицепленом файлике, на листе "Руками" - я простыми формулами посчитал - и там, в первом приближении - значение коэффициента Yt - в точности совпадает со значением из примера.

Цитата
panix1111 написал: Понял, чего точно не хватет - второй температуры, первая при которой меряли д ля расчета приближения, а вторая(для которой расчет и идет) должна использоваться на последующих шагах
Хм... Так вроде в алгоритме указано использовать именно замеренную температуру то есть PtP... Пошёл читать в сто первый раз алгоритм...

Вот же:
Цитата
Примечания    
При определении значения плотности p15, в формулы (1) и (3) подставляют значения температуры и давления нефти или нефтепродукта, при которых была измерена плотность PtP
НО
При расчете значения плотности PtP в формулы (1) и (3) подставляют значения температуры и давления нефти, при которых требуется определить плотность.
В рассматриваемом примере мы выводим именно p15

Цитата
Hypohelix написал: У меня возникло такое ощущение, что в стандарте где-то что-то пропущено. У меня в результате вычислений получаются те же цифры, что и у Вас. При этом, более-менее близкие к примеру цифры получаются только при убирании знака минус в формуле 1 у коэффициента Б15.  
Угу.... тоже игрался с этими минусами... И потом - вряд ли в госте напортачили... Я сверялся с предыдущими гостами - там практически те же самые формулы используются...
Изменено: Святослав И - 29.10.2017 17:25:11
 
хм, таки да
там в знаменателе - 0,99 (что немного увеличивает числитель), но можно принять за 1

а вот числитель действительно дает разницу (-) перед В15

(old_dens * Exp(b15 * (t - 15) * Abs((1 + 0.8 * b15 * (t - 15)))))
845.30955392501

(old_dens * Exp(-b15 * (t - 15) * Abs((1 + 0.8 * b15 * (t - 15)))))
827.08969661311

что то там не то с вычислениями
Код
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 - 29.10.2017 19:58:16
Мы в Екселе не работаем, мы в нём живём!
 
Цитата
panix1111 написал:
что то там не то с вычислениями
Угу, мягко говоря - не то... Настолько не то, что оно уходит в противоположном направлении. С другой стороны, вот ещё раз формула:


как вот её ещё можно понять иначе? У меня, грешным делом, была мысль, что скобки квадратные-фигурные каким-то образом могут влиять на очерёдность вычислений...
Изменено: Святослав И - 29.10.2017 20:07:13
 
Цитата
Святослав И написал:
скобки квадратные-фигурные каким-то образом могут влиять на очерёдность вычислений
квадратные - это модуль, я его в последнем варианте на всяк случай добавил в ф-лу (Abs), но всё равно что-то не сходится с расчетами из примера.
Что могу подсказать - есть готовые таблицы для расчета, может их как-то приспособить?
Мы в Екселе не работаем, мы в нём живём!
 
Цитата
panix1111 написал:
есть готовые таблицы для расчета, может их как-то приспособить?
Да, такой вариант, за неимением лучшего, я уже приспособил и что-то оно вычисляет. Но. Погрешности к сож., слишком велики. Велики на столько, что всю ценность вычислений сводят на "нет". Ну то есть, грубо говоря, на бензовозе в котором плещется 22-23 тонны солярки - литров 200 свободно "гуляет" туда/сюда только на этих знаках после запятой...  Соответственно, какие-то недоливы/сливы можно даже не пытаться ловить, соответственно - про какой-либо вменяемый учёт/контроль можно забывать... Просто я покурил форумы метрологов и там вполне чётко прослеживается мысль, что только вот эти формулы...
 
там есть такой момент:
Цитата
При определении значения плотности p15 в формулы (1) и (3) подставляют значения температуры и давления нефти или нефтепродукта, при которых была измерена плотность PtP. При расчете значения плотности PtP в формулы (1) и (3) подставляют значения температуры и давления нефти, при которых требуется определить плотность.
Изменено: panix1111 - 29.10.2017 21:22:43
Мы в Екселе не работаем, мы в нём живём!
 
panix1111, так я ровно о том же самом написал в #7-ом сообщении. Фишка в том, что мы как раз вычисляем именно р15, а не РtP... В алгоритме подстановок вроде об этом явно написано, или я как-то не так его читаю?

Цитата
3) Значение P15, вычисленное в первом приближении, подставляют в формулы (2) и (3) и вычисляют во втором приближении значения B15 и Yt, соответственно
Хотя вот начиная с этого пункта возникает некоторая неопределённость. Но опять же. Начиная с самой первой подстановки полученных b15 и Yt (значения которых, кстати, вычисляются чётко по примеру) - формула возвращает неверное значение, т.е. уже на самом первом шаге.
Изменено: Святослав И - 29.10.2017 21:56:30
 
Кросс:
http://www.programmersforum.ru/showthread.php?t=313704
http://www.excel-vba.ru/forum/index.php?topic=5355.0
http://forum.msexcel.ru/index.php/topic,11785.0.html
http://www.excelworld.ru/forum/10-35879-1
Я сам - дурнее всякого примера! ...
 
Святослав И, размещаете тему на нескольких форумах - информируйте об этом.
 
грубый хак, но расчет соответсвует гост-примеру:
Код
pTp = ((dens * Exp(-b15 * (t - 15) * Abs((1 + 0.8 * b15 * (t - 15))))) / (1 - Yt * p)) + 14.645

где 14.645 - просто ручками подобранный кофициэнт
Мы в Екселе не работаем, мы в нём живём!
 
Цитата
panix1111 написал:
14.645 - просто ручками подобранный кофициэнт
На других вводных, значения этого "коэффициента" будут другими. Что делает такое решение не вполне жизнеспособным...КМК... ))
 
Скрытый текст
Изменено: С.М. - 30.10.2017 21:43:01
 
С.М., А можно в двух словах буквально? Но в любом случае - спасибо. Опробую, отпишусь.
Изменено: Святослав И - 29.10.2017 23:03:53
 
Цитата
Святослав И написал:
А можно в двух словах буквально?
На одном листе:
 
С.М., см. лист "Анализ" - референсные данные из примера и рядом то, что выдают формулы.
Изменено: Святослав И - 30.10.2017 09:00:56
 
Нас нае... обманули. Расходимся. Оказывается (!!!) для того чтобы формула из стандарта, вот эта:


Работала правильно для приведения PtP в P15, её всего лишь (!!!) нужно перевернуть. Вот так:


И тогда всё внезапно начинает правильно работать. То, что в ГОСТе забыли об этом написать, ну так то мелочи. Нет слов.


Всем спасибо за посильный вклад!

P.S. Это мне народ на метрологическом форуме подсказал. Для них это - нечто само собой разумеющееся.
P.P.S. ...Как и здесь для отдельных людей - другие вещи.

Изменено: Святослав И - 01.11.2017 22:48:23
 
Святослав И, извините, а что Вы рассчитывали на листе Анализ ? Плотность при t=15 и P=0 или какую ?

Ещё, в РЕКОМЕНДАЦИИ написали:
"Значение плотности ρ15 находят методом последовательных приближений, используя итерационный метод прямых подстановок",
но забыли упомянуть, что ρ15, в этом случае определяется по "обратной" формуле:
Изменено: С.М. - 30.10.2017 12:07:15
 
Тьфу, картинка не появилась. Дублирую:
 
Цитата
С.М. написал:
но забыли упомянуть, что ρ15, в этом случае определяется по "обратной" формуле:
Я ничего не забыл. В стандарте об этом ни слова. К слову, я там выше привёл формулу, по которой оно сработало. Как и ваша.
 
Возвращаю файл (с вопросом):
 
Цитата
Святослав И написал #22:
Изменено: Святослав И  - 30 Окт 2017 12:10:18
Не, хорошо, что в реальном мире МашиныВремени - не существует !  :) .
 
Цитата
С.М. написал:
Возвращаю файл (с вопросом):
Т.е. ниже   Вы пытаетесь определить плотность опять при тех же 27.3°С, 2.45 МПа, и   измеренной плотности 836,15 кг/м3 ???
В файле, то что вы VBA проделываете с помощью цикла, я просто разложил на листе - для большей наглядности. Так лучше видно - на какой итерации какие результаты выходят. Собственно, debug.print в консоль (в вашей функции) приносит ровно те же результаты. То, что я пытался делать - привести плотность при 27,3°С в плотность при 15°С, при том же P (т.е. 2,45 МПа) для пересчёта в другое значение P - там есть свои формулы. При этом чётко следовал, рекомендациям, о том, что дескать нужные значения просто подставить в ту же формулу (1). То, что вычислять нужно на перевёрнутой формуле - ребята забыли упомянуть (равно как и привести её). Ну что ж, бывает...

Ещё раз, всем спасибо за посильное участие и помощь.
 
Цитата
Святослав И написал:
То, что вычислять нужно на перевёрнутой формуле - ребята забыли упомянуть
:)
https://www.youtube.com/watch?v=Q-KCknBxVtg
 
Аааааа ! Нашёл у себя в коде ошибку - запутался в скобках (и на старухъ бывает прорух).
Заменил в #18, и новая редакция файла:
Изменено: С.М. - 30.10.2017 21:42:01
Страницы: 1 2 След.
Наверх