solving a system of equations where some variables are linked to a function

I'm analysing a thermodynamic model where I must solve a system of equations for a certain component.

The enthalpy however depends on the temperature and the concentration, which are also unkown variables. is it possible to add the function of the relation between the enthalpy and the concentration in the system of equations?

See below the script that I tried. (eqn 7 and 8 are the ones that don't work currently)

%%variables that are known.

m_3=0.004 ;

h_3=1530;

y_3=0.62;

To=25+273;

P_c=1.5 ;

UA=0.5;

hfg2=2256;

hfg=1815;

%%determine unknown variables.

syms x_5 y_6 m_5 m_6 h_5 h_6 Q_rect T

eqn1= m_3*h_3-m_5*h_5-m_6*h_6-Q_rect==0;

eqn2=Q_rect-m_5*(x_5*hfg+(1-x_5)*hfg2)==0;

eqn3=Q_rect-UA*(T-To)==0;

eqn4=m_3*y_3-m_5*x_5-m_6*y_6==0;

eqn5=m_3-m_6-m_5==0;

eqn6=m_3*(1-y_3)-(1-x_5)*m_5==0;

eqn7=h_5==enthalpyliquid(x_5,T);

eqn8=h_6==enthalpygas(y_6,P_c);

[solm5 ,solm6, solh5, solh6 ,Qrect ,T]=vpasolve([eqn1,eqn2,eqn3,eqn4,eqn5,eqn6,eqn7,eqn8],[x_5, y_6,m_5,m_6,h_5,h_6,Q_rect,T,h_5ex,h_6ex]);

[darova]

Any chance to see these?

[Emma Van Puyenbroeck]

function hl= enthalpyliquid (y,T)

m=[0 0 0 0 0 0 1 1 2 3 5 5 5 6 6 8];

n=[1 4 8 9 12 14 0 1 1 3 3 4 5 2 4 0];

a= [-7.61080 25.6905 -247.092 325.952 -158.854 61.9084 11.4314 1.18157 2.84179 7.41609 891.844 -1613.09 622.106 -207.588 -6.87393 3.50716];

h0=100 ;

TO=273.16;

h2=0;

for k=1:16

h1=a(k).*((T./TO)-1).^m(k).*y.^n(k);

h2=h2+h1;

end

hl=h0*h2;

[Emma Van Puyenbroeck]

function hg=enthalpygas(x,T)

m=[0 1 2 3 0 1 2 3 0 1 2 0 1 0 4 2 1];

n=[0 0 0 0 2 2 2 2 3 3 3 4 4 5 6 7 10];

a=[1.28827 0.125247 -2.08748 2.17696 2.35687 -8.86987 10.2635 -2.37440 -6.70515 16.4508 -9.36849 8.42254 -8.58807 -2.77049 -0.961248 0.988009 0.308482];

h0=1000;

TO=324;

h2=0;

for k=1:17

h1= a(k).*(1-(T/TO)).^m(k).*(1-x).^(n(k)/4);

h2=h1+h2;

end

hg=h0*h2;

[Emma Van Puyenbroeck]

Would you maybe recommend another solver?

[Walter Roberson]

Those two functions do not do any logical tests on the inputs, and do not use mod() or numeric integration or any complicated functions, and only calculate deterministic values. Under those circumstances, you can pass symbolic variables to the functions to calculate a single formula that expresses the output of the function in terms of the two inputs. The effect is that were you call the function with symbolic inputs it is as-if you replaced the call with the appropriate long formula in construction of the equations, replacing the function call with the expression the function calculates. At that point it stops mattering that you have written a function call: you could have written the messy formula instead right in the equation.

So as far as MATLAB is concerned, under those stated conditions, there is no difference between whether you call a function or write the expression. There is no difference in whether MATLAB can solve the system or not. You really just have a set of simultaneous equations, possibly nonlinear

[Emma Van Puyenbroeck]

But I receive the error

'Error using mupadengine/feval (line 187)

Symbolic parameters not supported in nonpolynomial equations. '

So how should I implement this expression or function? Or should I use a solver that can work with nonlinear equations?

Answer 1

[darova]

I tried fsolve

%%variables that are known. 
m_3=0.004 ;
h_3=1530;
y_3=0.62; 
To=25+273; 
P_c=1.5 ; 
UA=0.5;
hfg2=2256; 
hfg=1815;
% x_5  y_6  m_5  m_6  h_5  h_6  Q_rect  T 
% x(1) x(2) x(3) x(4) x(5) x(6) x(7)    x(8)
F = @(x) [m_3*h_3-x(3)*x(5)-x(4)*x(6)-x(7)
            x(7)-x(3)*(x(1)*hfg+(1-x(1))*hfg2)
            x(7)-UA*(x(8)-To)
            m_3*y_3-x(3)*x(1)-x(4)*x(2)
            m_3-x(4)-x(3)
            m_3*(1-y_3)-(1-x(1))*x(3)
            x(5)-enthalpyliquid(x(1),x(8))
            x(6)-enthalpygas(x(2),P_c)];
opt1 = optimoptions('fsolve',...
                    'Display','iter',...
                    'TolFun',1e-6,...
                    'TolX',1e-6,...
                    'MaxIter',1000,...
                    'MaxFunEvals',6500);

[X,fval] = fsolve(F,ones(1,8)*1,opt1);
X'
fval