Основываясь на ваших комментариях, я собрал быстрый сценарий, который поможет вам начать работу в правильном направлении - в дополнение к нескольким заметкам.
Во-первых, то, что вы пытаетесь решить, является нетривиальным, особенно учитывая, что Cp (T) сильно нелинейна и сверху этого расходится вблизи критической точки.
Кроме того, учитывая, что H = CpdT, и вы используете REFPROP, более точным вычислением будет H (T_out) - H (T_in), который будет включать температурные эффекты.
С учетом сказанного, выписав свои балансы массы, вы можете сформулировать систему нелинейных уравнений для решения, как показано ниже. Поскольку я не мог решить ваше использование UA
, я просто заменил его на «коэффициент эффективности». Не стесняйтесь использовать свое соглашение о выборе (например, NTU), чтобы заполнить его.
Я предупреждаю вас, что то, что я разместил ниже, является неполным. Это отправная точка и есть несколько вопросов, которые необходимо учитывать:
Я с указанием n-2
переменных температурами, обеспечивая при этом n
уравнения, означающем, что границы могут быть выполнено в виде Ta_end = Tb_end (вы увидите, что я имею в виду, если вы бежите это и сюжет). Это также может привести к разрыву температуры в зависимости от стадии.
В зависимости от ваших давлений и жидкостей, если вы проходите фазовый переход 1-го порядка, вы можете столкнуться с некоторыми странными численными проблемами.
Если вы можете принять постоянную теплоемкость, вы можете переформулировать матрицу ниже в гораздо более легкую и разрешимую систему. Я предлагаю этот более сложный пример в качестве отправной точки.
function Tf = gen_mat
n = 10; % Number of stages
A = zeros(n,2*n+2); % Coefficients matrix.
ma = 10; % Mass flow of A.
mb = 5; % Mass flow of B.
eff = 1.0; % Effectiveness factor.
Te(1) = 383.15; % Inlet T of A.
Te(2) = 278.15; % Inlet T of B.
P = 15000; % P (kPa)
% Generate guesses
Ti = linspace(Te(1),Te(2),n+1)';
Tg = zeros(2*n+1,1);
for i=1:n+1
Tg(2*i-1) = Ti(i);
Tg(2*i) = Ti(i);
end
Tg([1,2*n+1]) = []; % We are only interested in the middle sections.
% construct coefficient matrix.
for i=1:n
A(i,2*i-1) = -ma;
A(i,2*i) = eff*mb;
A(i,2*i+1) = ma;
A(i,2*i+2) = -eff*mb;
end
% solve system of nonlinear equations
Tf = fsolve(@(T)obj_fxn(T,A,P,Te), Tg, optimset('Display', 'iter'));
Tf = [Te(1);Tf;Te(2)]; % Append temperatures.
end
function b = obj_fxn(T,A,P,Te)
T = [Te(1); T; Te(2)];
for i=1:2:length(T)-1
x(i) = refpropm('H','T',T(i),'P',P,'CO2');
x(i+1) = refpropm('H','T',T(i+1),'P',P,'water');
end
b = A*x';
end
Вы можете создать сырой участок выше после запуска это следующим образом: plot(0:10,Tf(1:2:end)-273,10:-1:0,Tf(2:2:end)-273)
.
Извините, что это не более организованный, но я надеюсь, что это поможет.
Сохраняется ли постоянная константа 'T_a_in' и' T_b_in' ** ** или они изменяются, когда вы пытаетесь найти решение для системы? – rayryeng
Как определяются функции 'cp'? Это, вероятно, будет иметь важное значение. – Jubobs
В дополнение к тому, что спросил Джубос, если бы вы могли также предоставить точные определения ваших функций 'cp' (код и все), это было бы полезно для понимания и, возможно, решения вашей проблемы. – rayryeng