1. Обязательно представиться на русском языке кириллицей (заполнить поле "Имя").
  2. Фиктивные имена мы не приветствуем. Ивановых и Пупкиных здесь уже достаточно.
  3. Не писать свой вопрос в первую попавшуюся тему - вместо этого создать новую тему.
  4. За поиск, предложение и обсуждение пиратского ПО и средств взлома - бан без предупреждения.
  5. Рекламу и частные объявления "куплю/продам/есть халтура" мы не размещаем ни на каких условиях.
  6. Перед тем как что-то написать - читать здесь и здесь, а студентам - обязательно здесь.
  7. Не надо писать в ЛС администраторам свои технические вопросы. Администраторы форума отлично знают как работает форум, а не все-все контроллеры, о которых тут пишут.

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

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

Автор темы
VADR
администратор
администратор
Сообщения: 4711
Зарегистрирован: 25 июл 2008, 07:12
Имя: Диев Александр Васильевич
Страна: Россия
город/регион: г. Сегежа, Карелия
Благодарил (а): 192 раза
Поблагодарили: 336 раз

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

Сообщение VADR »

Недавно довелось столкнуться с такой задачей - расчёт тепловой энергии пара и воды в разных диапазонах давлений и температур. Помня такую книгу, как "Термодинамические свойства воды и водяного пара", Вукалович М.П., предполагал, что тут всё непросто и аппроксимировать кривые можно до потери пульса. Но откуда-то из закоулков памяти вылезло и другое: существуют приближённые (тем не менее, весьма точные) формулы расчёта. По информации, которую дал гугль, разработал эти формулы "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
почётный участник форума
почётный участник форума
Сообщения: 1746
Зарегистрирован: 18 янв 2009, 12:25
Имя: Дмитрий
Страна: Россия
город/регион: Москва
Благодарил (а): 6 раз
Поблагодарили: 132 раза

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

Сообщение hell_boy »

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

Автор темы
VADR
администратор
администратор
Сообщения: 4711
Зарегистрирован: 25 июл 2008, 07:12
Имя: Диев Александр Васильевич
Страна: Россия
город/регион: г. Сегежа, Карелия
Благодарил (а): 192 раза
Поблагодарили: 336 раз

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

Сообщение VADR »

hell_boy писал(а): существуют также ГСССД 187-89 и МИ 2451-98
Блин, не знал.
hell_boy писал(а): Макросы VB лежали на сайте IAPWS
А этого дела что-то так и не нашёл. Или плохо искал. Знаю, что есть софтина под названием WaterSteamPro (и у них даже имена функций начинаются так же, как у меня - с wsp, однако под последней "p" я имел в виду "properties", а не "pro") и она сделана на основе того же самого документа, но в моём случае применить её вряд ли было возможно.
Повторное использование кода не отменяет повторного использования мозга при его повторном использовании.
Ответить

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