система нелинейных уравнений с экспоненциальным матлабом - PullRequest
Новая
       64

система нелинейных уравнений с экспоненциальным матлабом

0 голосов
/ 27 мая 2020

Я пытаюсь решить систему из двух нелинейных уравнений, одно экспоненциальное, с двумя неизвестными. Я пробовал использовать fsolve, но система очень нестабильна, и я не понимаю, как заставить ее работать. Я понял, что если мощность равна нулю, алгоритм работает, но с любым другим числом результаты сильно колеблются. CoolProp - это библиотека, которую можно найти в Интернете.

Не могли бы вы рассказать мне, почему это произошло, или не могли бы вы предложить мне более подходящую функцию для решения проблемы?

Большое спасибо 

clc;clear;
%data H2
Ru=8.314; 
MM_H2=2.016*1e-3; 
Rgas=Ru/MM_H2; 
beta=1.9155*1e-6;
V_vt=1.288; 
T0=288; 
p0_vt=20*1e-1; 
T8=288; 
Tf=288; 
m0=2.14; 
A_in=1.5;
h_conv_in=200; 
APRR=0.003;

t=1:1:60000;

cp_vt=zeros(1,length(t))';
cv_vt=zeros(1,length(t))';
gamma_vt=zeros(1,length(t))';
m_vt=zeros(1,length(t))';
m_flow=zeros(1,length(t))';
T_vt=zeros(1,length(t))';
P_vt=zeros(1,length(t))';



for i=1:length(t)

    %to calculate hydrogen pressure during refueling
    P_vt(i)=p0_vt+APRR*t(i); %[MPa]

    if i==1
        %to calculate hydrogen mass and temperature in vt
        cp_vt(i)=CoolProp.PropsSI('C', 'P', p0_vt*1e6, 'T', T0, 'H2');
        cv_vt(i)=CoolProp.PropsSI('O', 'P', p0_vt*1e6, 'T', T0, 'H2');
        gamma_vt(i)=cp_vt(i)/cv_vt(i); %heat capacity ratio

        %write the two functions
        F = @(x) [P_vt(i) * V_vt - x(1) * Rgas*1e-6*x(2)*(1+beta*1e6*P_vt(i)/x(2));
            x(2)-(m0/x(1))^(1+h_conv_in*A_in/((x(1)-m0)*cv_vt(i)))*T0-...
            (1-(m0/x(1))^(1+h_conv_in*A_in/((x(1)-m0)*cv_vt(i))))*(gamma_vt(i)*T8+Tf*h_conv_in*A_in/((x(1)-m0)*cv_vt(i)))/...
            (1+h_conv_in*A_in/((x(1)-m0)*cv_vt(i)))];

        x0 = [2.15;288];
        options=optimset('Display','off','TolX', 1e-4);

        [x,fval] = fsolve(F,x0,options);
        m_vt(i)=x(1);
        T_vt(i)=x(2);
        m_flow(i)=m_vt(i)-m0;

    else
        %to calculate hydrogen mass and temperature in vt
        cp_vt(i)=CoolProp.PropsSI('C', 'P', P_vt(i-1)*1e6, 'T', T_vt(i-1), 'H2');
        cv_vt(i)=CoolProp.PropsSI('O', 'P', P_vt(i-1)*1e6, 'T', T_vt(i-1), 'H2');
        gamma_vt(i)=cp_vt(i)/cv_vt(i); %heat capacity ratio


        x0 = [5;290];
        options=optimset('Display','off','TolX', 1e-4);

        F = @(x) [P_vt(i) * V_vt - x(1) * Rgas*1e-6*x(2)*(1+beta*1e6*P_vt(i)/x(2));
            x(2)-(m0/x(1))^(1+h_conv_in*A_in/((x(1)-m_vt(i-1))*cv_vt(i)))*T0-...
            (1-(m0/x(1))^(1+h_conv_in*A_in/((x(1)-m_vt(i-1))*cv_vt(i))))*(gamma_vt(i)*T8+Tf*h_conv_in*A_in/((x(1)-m_vt(i-1))*cv_vt(i)))/...
            (1+h_conv_in*A_in/((x(1)-m_vt(i-1))*cv_vt(i)))];

        [x,fval] = fsolve(F,x0,options);
        m_vt(i)=x(1);
        T_vt(i)=x(2);
        m_flow(i)=m_vt(i)-m_vt(i-1);
    end
end
...