Профиль должен быть заполнен на русском языке кириллицей. Заполнение профиля заведомо ложными или некорректными данными - причина возможного отказа в регистрации на форуме.

Расчёт свойств воды и водяного пара

Аватара пользователя

Автор темы
VADR
администратор
администратор
Сообщения: 2411
Зарегистрирован: 25 июл 2008, 06:12
Ф.И.О.: Диев Александр Васильевич
Благодарил (а): 19 раз
Поблагодарили: 26 раз

Расчёт свойств воды и водяного пара

Сообщение VADR » 12 мар 2016, 01:01

Недавно довелось столкнуться с такой задачей - расчёт тепловой энергии пара и воды в разных диапазонах давлений и температур. Помня такую книгу, как "Термодинамические свойства воды и водяного пара", Вукалович М.П., предполагал, что тут всё непросто и аппроксимировать кривые можно до потери пульса. Но откуда-то из закоулков памяти вылезло и другое: существуют приближённые (тем не менее, весьма точные) формулы расчёта. По информации, которую дал гугль, разработал эти формулы "The International Association for the Properties of Water and Steam", текст существует в открытом доступе и лежит здесь: http://www.iapws.org/relguide/IF97-Rev.pdf. Документ, по сравнению с Вукаловичем, "всего" на 49 страниц :).
Формулы на первый взгляд кажутся многоэтажными, но это только на первый взгляд. Важно точно вбить таблицы (и для этого язык, на котором вы будете это реализовывать, должен поддерживать массивы), а сами формулы... впрочем, ниже я приведу то, что у меня получилось - видно, что ничего принципиально сложного нет.
В моей ситуации я был ограничен в выборе языка программирования: целевая система поддерживает расчёты:
1. На своём внутреннем скриптовом языке (не поддерживает массивов и циклов).
2. В виде формул в Excel (назначаются ячейки для входящих данных и результата, можно прописать любые расчёты, но тупит это просто неимоверно - раз в секунду вызывать не получится никак).
3. COM-объект. Даже не стал влезать в это дело: мало того, что саму программу надо написать, её ещё надо специальным образом зарегистрировать, о ней и её регистрации надо не забывать при резервном копировании... Я на такие жертвы пойду только в том случае, когда других вариантов не останется совсем.
4. EXE-файл. Что-то типа предыдущего пункта, так что тоже - самый крайний случай.
5. VBScript. Массивы и циклы поддерживаются, текст функций хранится непосредственно в базе данных целевой системы и бэкапится вместе с ней без каких-либо дополнительных действий, текст функций всегда виден во встроенном редакторе.
В общем, хоть я и недолюбливаю Visual Basic вообще и VBScript в частности, в данном случае пришлось остановиться именно на нём. Это я тут налепил лирическое отступление, чтобы потом не отвечать на вопросы, почему именно VBS.
В общем, переходя от лирики к технике:
1. Из всех рассчитываемых параметров мне нужны были только удельный объём и энтальпия для воды и пара. Если посмотреть на функции в первоисточнике, нет ничего сложного на базе этих функций сделать аналогичные для остальных параметров (энтропия, внутренняя энергия, изохорная и изобарная теплоёмкости, скорость звука)
2. Для линии насыщения - нужны зависимость температуры от давления и наоборот.
3. Если есть возможность, следует использовать переменные двойной точности (double). По крайней мере, в тестовых примерах так рекомендуется.
4. На странице 4 документа есть график, где отмечены зоны, в которых применима та или иная формула. В моём случае в технологии нет пара и воды с давлением выше ~16МПа либо температурой выше 1073,15K, так что меня во всех этой красоте интересуют только регионы 1 и 2 и линия насыщения 4. Принцип расчёта в других регионах ничем не отличается, так что желающие могут реализовать эти формулы самостоятельно.
Итак, что у меня получилось:
1. Удельный объём воды (там, где температура ниже линии насыщения), по графику - регион 1:

Код: Выделить всё

Function wspV1(ByVal P, ByVal T)
   Dim Pi, Tau, RES, n, I, J, k, R
   R = 0.461526
   n = Array (0, 0.14632971213167, -0.84548187169114, -0.37563603672040E1, 0.33855169168385E1, -0.95791963387872, 0.15772038513228, -0.16616417199501E-1, 0.81214629983568E-3, 0.28319080123804E-3, -0.60706301565874E-3, -0.18990068218419E-1, -0.32529748770505E-1, -0.21841717175414E-1, -0.52838357969930E-4, -0.47184321073267E-3, -0.30001780793026E-3, 0.47661393906987E-4, -0.44141845330846E-5, -0.72694996297594E-15, -0.31679644845054E-4, -0.28270797985312E-5, -0.85205128120103E-9, -0.22425281908000E-5, -0.65171222895601E-6, -0.14341729937924E-12, -0.40516996860117E-6, -0.12734301741641E-8, -0.17424871230634E-9, -0.68762131295531E-18, 0.14478307828521E-19, 0.26335781662795E-22, -0.11947622640071E-22, 0.18228094581404E-23, -0.93537087292458E-25)
   I = Array (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32)
   J = Array (0, -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41)
    Pi = P / 16.53
   Tau = 1386 / T
   RES = 0
   for k = 1 to 34
      RES = RES + (-1 * n(k) * I(k) * ((7.1 - Pi) ^ (I(k)-1)) * ((Tau - 1.222)^J(k)))
   Next
   wspV1 = RES * Pi * R * T / (P * 1000)
end Function

2. Удельная энтальпия воды:

Код: Выделить всё

Function wspH1(ByVal P, ByVal T)
   Dim Pi, Tau, RES, n, I, J, k, R
   R = 0.461526
   n = Array (0, 0.14632971213167, -0.84548187169114, -0.37563603672040E1, 0.33855169168385E1, -0.95791963387872, 0.15772038513228, -0.16616417199501E-1, 0.81214629983568E-3, 0.28319080123804E-3, -0.60706301565874E-3, -0.18990068218419E-1, -0.32529748770505E-1, -0.21841717175414E-1, -0.52838357969930E-4, -0.47184321073267E-3, -0.30001780793026E-3, 0.47661393906987E-4, -0.44141845330846E-5, -0.72694996297594E-15, -0.31679644845054E-4, -0.28270797985312E-5, -0.85205128120103E-9, -0.22425281908000E-5, -0.65171222895601E-6, -0.14341729937924E-12, -0.40516996860117E-6, -0.12734301741641E-8, -0.17424871230634E-9, -0.68762131295531E-18, 0.14478307828521E-19, 0.26335781662795E-22, -0.11947622640071E-22, 0.18228094581404E-23, -0.93537087292458E-25)
   I = Array (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32)
   J = Array (0, -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41)
    Pi = P / 16.53
   Tau = 1386 / T
   RES = 0
   for k = 1 to 34
      RES = RES + ( n(k) * ((7.1 - Pi) ^ I(k)) * J(k) * ((Tau - 1.222)^(J(k)-1)))
   Next
   wspH1 = RES * Tau * R * T
end Function

3. Удельный объём пара (температура выше линии насыщения), по графику - регион 2:

Код: Выделить всё

Function wspV2(ByVal P, ByVal T)
   Dim n0, J0, RES, R, Pi, Tau, k, n, I, J
   R = 0.461526
   n0 = Array (0, -0.96927686500217E1, 0.10086655968018E2, -0.56087911283020E-2, 0.71452738081455E-1, -0.40710498223928, 0.14240819171444E1, -0.43839511319450E1, -0.28408632460772, 0.21268463753307E-1)
    J0 = Array (0, 0, 1, -5, -4, -3, -2, -1, 2, 3)
   I = Array (0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24)
   J = Array (0, 0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58)
   n = Array (0, -0.17731742473213E-2, -0.17834862292358E-1, -0.45996013696365E-1, -0.57581259083432E-1, -0.50325278727930E-1, -0.33032641670203E-4, -0.18948987516315E-3, -0.39392777243355E-2, -0.43797295650573E-1, -0.26674547914087E-4, 0.20481737692309E-7, 0.43870667284435E-6, -0.32277677238570E-4, -0.15033924542148E-2, -0.40668253562649E-1, -0.78847309559367E-9, 0.12790717852285E-7, 0.48225372718507E-6, 0.22922076337661E-5, -0.16714766451061E-10, -0.21171472321355E-2, -0.23895741934104E2, -0.59059564324270E-17, -0.12621808899101E-5, -0.38946842435739E-1, 0.11256211360459E-10, -0.82311340897998E1, 0.19809712802088E-7, 0.10406965210174E-18, -0.10234747095929E-12, -0.10018179379511E-8, -0.80882908646985E-10, 0.10693031879409, -0.33662250574171, 0.89185845355421E-24, 0.30629316876232E-12, -0.42002467698208E-5, -0.59056029685639E-25, 0.37826947613457E-5, -0.12768608934681E-14, 0.73087610595061E-28, 0.55414715350778E-16, -0.94369707241210E-6)
   Pi = P
   Tau = 540 / T
   RES = 1 / Pi
   For k = 1 to 43
      RES = RES + (n(k) * I(k) * (Pi ^ (I(k) - 1)) * ((Tau - 0.5) ^ J(k)))
   Next
   wspV2 = RES * Pi * R * T / (P * 1000)
end Function

4. Удельная энтальпия пара:

Код: Выделить всё

Function wspH2(ByVal P, ByVal T)
   Dim n0, J0, RES, R, Pi, Tau, k, n, I, J
   R = 0.461526
   n0 = Array (0, -0.96927686500217E1, 0.10086655968018E2, -0.56087911283020E-2, 0.71452738081455E-1, -0.40710498223928, 0.14240819171444E1, -0.43839511319450E1, -0.28408632460772, 0.21268463753307E-1)
    J0 = Array (0, 0, 1, -5, -4, -3, -2, -1, 2, 3)
   I = Array (0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 10, 10, 10, 16, 16, 18, 20, 20, 20, 21, 22, 23, 24, 24, 24)
   J = Array (0, 0, 1, 2, 3, 6, 1, 2, 4, 7, 36, 0, 1, 3, 6, 35, 1, 2, 3, 7, 3, 16, 35, 0, 11, 25, 8, 36, 13, 4, 10, 14, 29, 50, 57, 20, 35, 48, 21, 53, 39, 26, 40, 58)
   n = Array (0, -0.17731742473213E-2, -0.17834862292358E-1, -0.45996013696365E-1, -0.57581259083432E-1, -0.50325278727930E-1, -0.33032641670203E-4, -0.18948987516315E-3, -0.39392777243355E-2, -0.43797295650573E-1, -0.26674547914087E-4, 0.20481737692309E-7, 0.43870667284435E-6, -0.32277677238570E-4, -0.15033924542148E-2, -0.40668253562649E-1, -0.78847309559367E-9, 0.12790717852285E-7, 0.48225372718507E-6, 0.22922076337661E-5, -0.16714766451061E-10, -0.21171472321355E-2, -0.23895741934104E2, -0.59059564324270E-17, -0.12621808899101E-5, -0.38946842435739E-1, 0.11256211360459E-10, -0.82311340897998E1, 0.19809712802088E-7, 0.10406965210174E-18, -0.10234747095929E-12, -0.10018179379511E-8, -0.80882908646985E-10, 0.10693031879409, -0.33662250574171, 0.89185845355421E-24, 0.30629316876232E-12, -0.42002467698208E-5, -0.59056029685639E-25, 0.37826947613457E-5, -0.12768608934681E-14, 0.73087610595061E-28, 0.55414715350778E-16, -0.94369707241210E-6)
   Pi = P
   Tau = 540 / T
   RES = 0
   For k = 1 to 9
      RES = RES + n0(k) * J0(k) * (Tau ^ (J0(k) - 1))
   Next
   For k = 1 to 43
      RES = RES + (n(k) * (Pi ^ I(k)) * J(k) * ((Tau - 0.5) ^ (J(k) - 1)))
   Next
   wspH2 = RES * Tau * R * T
end Function

5. Температура на линии насыщения в зависимости от давления и наоборот:

Код: Выделить всё

Function wspTsat(ByVal P)
   Dim Beta, D, E, F, G, n
   n = Array (0, 0.11670521452767E4, -0.72421316703206E6, -0.17073846940092E2, 0.12020824702470E5, -0.32325550322333E7, 0.14915108613530E2, -0.48232657361591E4, 0.40511340542057E6, -0.23855557567849, 0.65017534844798E3)
   Beta = Sqr( Sqr (P))
   E = Beta^2 + (n(3)*Beta) + n(6)
   F = (n(1)*Beta^2) + (n(4)*Beta) + n(7)
   G = (n(2)*Beta^2) + (n(5)*Beta) + n(8)
   D = 2*G / (-1*F - Sqr(F^2 - 4*E*G))
    wspTsat = (n(10) + D - Sqr((n(10)+D)^2 - 4*(n(9)+n(10)*D)))/2
end Function

Код: Выделить всё

Function wspPsat(ByVal T)
   Dim Theta, A, B, C, n
   n = Array (0, 0.11670521452767E4, -0.72421316703206E6, -0.17073846940092E2, 0.12020824702470E5, -0.32325550322333E7, 0.14915108613530E2, -0.48232657361591E4, 0.40511340542057E6, -0.23855557567849, 0.65017534844798E3)
   Theta = T + (n(9) / (T - n(10)))
   A = Theta^2 + (n(1) * Theta) + n(2)
   B = (n(3) * Theta^2) + (n(4) * Theta) + n(5)
   C = (n(6) * Theta^2) + (n(7) * Theta) + n(8)
   wspPsat = (2*C/(-1*B + Sqr(B^2 - 4*A*C))) ^ 4
end Function


Единицы измерения - можно в первоисточнике посмотреть. Там же есть тестовые данные для расчёта. Кстати, в формулах для удельного объёма в предпоследней строке я дополнительно разделил результат на 1000, хотя в описании формул такого не было. В общем, есть там какая-то нестыковка с единицами измерения, в итоге, чтобы совпали результаты по тестовым выражениям, пришлось их поделить на 1000 :).

Вот как-то так, собственно. Кому вдруг понадобится - пользуйтесь на здоровье.
Повторное использование кода не отменяет повторного использования мозга при его повторном использовании.

Аватара пользователя

hell_boy
почётный участник форума
почётный участник форума
Сообщения: 1078
Зарегистрирован: 18 янв 2009, 12:25
Ф.И.О.: Дмитрий
Благодарил (а): 2 раза
Поблагодарили: 22 раза

Расчёт свойств воды и водяного пара

Сообщение hell_boy » 12 мар 2016, 14:09

Кроме IAPWS IF-97 существуют также ГСССД 187-89 и МИ 2451-98. Макросы VB лежали на сайте IAPWS
steam2010.xls
У вас нет необходимых прав для просмотра вложений в этом сообщении.
"Умные люди обсуждают идеи, средние - события, а глупые - людей" Л.Н. Толстой

Аватара пользователя

Автор темы
VADR
администратор
администратор
Сообщения: 2411
Зарегистрирован: 25 июл 2008, 06:12
Ф.И.О.: Диев Александр Васильевич
Благодарил (а): 19 раз
Поблагодарили: 26 раз

Расчёт свойств воды и водяного пара

Сообщение VADR » 12 мар 2016, 23:14

hell_boy писал(а):Источник цитаты существуют также ГСССД 187-89 и МИ 2451-98

Блин, не знал.
hell_boy писал(а):Источник цитаты Макросы VB лежали на сайте IAPWS

А этого дела что-то так и не нашёл. Или плохо искал. Знаю, что есть софтина под названием WaterSteamPro (и у них даже имена функций начинаются так же, как у меня - с wsp, однако под последней "p" я имел в виду "properties", а не "pro") и она сделана на основе того же самого документа, но в моём случае применить её вряд ли было возможно.
Повторное использование кода не отменяет повторного использования мозга при его повторном использовании.


Вернуться в «Заготовки для базы знаний»



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 0 гостей