2015-03-16 4 views
0

У меня есть система уравнений, определяющая систему. Я не знаю, как это сделать в Matlab, кроме того, что это возможно, и требуются итерации.Как решить набор уравнений Matlab

У меня есть три уравнения с тремя неизвестными (Q, T_a_out и T_b_out).

Система разделена на несколько частей. Информация о температуре в первой и последней подчасти должна использоваться для нахождения температур в подчасти в середине. Мне нужны Q, T_a_out и T_b_out для каждой подчасти.

System description

Уравнения, используемые для описания системы заключаются в следующем:

Q=U*((T_a_out-T_b_in)+(T_a_in-T_b_out))/2 
Q=m_a*(cp_a_in*T_a_in-cp_a_out*T_a_out) 
Q=m_b*cp_b*(T_b_out-T_b_in) 

Эти параметры известны:

Initial T_a_in (110) 
Initial T_b_in (5) 
U 
m_a 
m_b 
n (number of sub-parts) 
cp_b 

Значение cp является функцией температуры она принадлежит до:

cp_a_in is a function of the temperature: cp_a_in=function_a(T_a_in) 
cp_a_out is a function of the temperature: cp_a_out=function_a(T_a_out) 

Начальное значение a (110) имеет более высокое значение, чем конечная температура b_out в последней подчасти. Начальное значение b_in (5) имеет меньшее значение, чем конечное значение a_out.

Как рассчитать температуру out для каждой детали в Matlab?

+0

Сохраняется ли постоянная константа 'T_a_in' и' T_b_in' ** ** или они изменяются, когда вы пытаетесь найти решение для системы? – rayryeng

+1

Как определяются функции 'cp'? Это, вероятно, будет иметь важное значение. – Jubobs

+0

В дополнение к тому, что спросил Джубос, если бы вы могли также предоставить точные определения ваших функций 'cp' (код и все), это было бы полезно для понимания и, возможно, решения вашей проблемы. – rayryeng

ответ

0

Основываясь на ваших комментариях, я собрал быстрый сценарий, который поможет вам начать работу в правильном направлении - в дополнение к нескольким заметкам.

Во-первых, то, что вы пытаетесь решить, является нетривиальным, особенно учитывая, что 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).

Извините, что это не более организованный, но я надеюсь, что это поможет.

Смежные вопросы