Water Thermodynamic: pressure calculation

3 views (last 30 days)
Nico
Nico on 3 Dec 2021
Hi everyone, I'm writing the code of the water's pressure.
I wrote the matrix of the constant:
A= [ A11 A12 A13 A14 A15 A16 A17;
A21 A22 A23 A24 A25 A26 A27;
A31 A32 A33 A34 A35 A36 A37;
A41 A42 A43 A44 A45 A46 A47;
A51 A52 A53 A54 A55 A56 A57;
A61 A62 A63 A64 A65 A66 A67;
A71 A72 A73 A74 A75 A76 A77;
A81 A82 A83 A84 A85 A86 A87;
A91 A92 A93 A94 A95 A96 A97;
A101 A102 A103 A104 A104 A106 A107; ];
I want to write the same C++ code in Matlab:
double water::C(int i)
{
double tau = Ta/T;
return (i == 0 ? R*T : R*T*(tau - tauc)*pow(tau - taua[i],i-1));
}
double water::H(int j)
{
double factor, sum, rho_aj;
rho_aj = (j == 0 ? Roa1 : Roaj);
sum = 0.0;
factor = Rho - rho_aj;
for (int i=6; i>0; i--) {
sum += (A[i][j] + Rho*(i+1)*A[i+1][j]);
sum *= factor;
}
sum += (A[0][j] + Rho*A[1][j]);
sum += (exp(-E*Rho)*((1.0 - Rho*E)*A[8][j]
+ Rho*(2.0 - Rho*E)*A[9][j]));
sum += A[7][j]*pow(factor,7);
return Rho*Rho*sum;
}
double water::Pp()
{
double P = Rho*R*T;
for (int i=0; i<7; i++) {
P += C(i)*H(i);
}
return P;
}
These are the 3 function: the last one is the pressure.
I wrote the H function like this:
function [H] = H_fun(rho, j)
global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
u0 s0 M Tc pc rhoC T0 tauC
A = [ A11 A12 A13 A14 A15 A16 A17;
A21 A22 A23 A24 A25 A26 A27;
A31 A32 A33 A34 A35 A36 A37;
A41 A42 A43 A44 A45 A46 A47;
A51 A52 A53 A54 A55 A56 A57;
A61 A62 A63 A64 A65 A66 A67;
A71 A72 A73 A74 A75 A76 A77;
A81 A82 A83 A84 A85 A86 A87;
A91 A92 A93 A94 A95 A96 A97;
A101 A102 A103 A104 A104 A106 A107; ];
if j==1
rho_Aj = rhoA1;
else
rho_Aj = rhoAj;
end
sum = 0.0;
factor = rho - rho_Aj;
for i = [ 7 6 5 4 3 2 1]
sum = sum + (A(i,j) + rho*(i+1)*A((i+1),j));
sum = sum*factor;
end
sum = sum + (A(1,j) + rho*A(2,j));
sum = sum + ( exp(-E*rho)*((1.0 - rho*E)*A(9,j) + rho*(2.0 - rho*E)*A(10,j)) );
sum = sum + A(8,j)*(factor^(7));
H = rho*rho*sum;
end
The C function like this:
function [C]= C_fun(T, j)
global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
u0 s0 M Tc pc rhoC T0 tauC
tauAj = [ 1.544912 2.5 2.5 2.5 2.5 2.5 2.5 ];
tau = Ta / T;
if j==1
C = R*T;
else
C = R*T * (tau - tauC)*((tau - tauAj(j))^(j-2));
end
and the Pressure function:
function [P] = P_fun(rho, T)
global R A11 A21 A31 A41 A51 A61 A71 A81 A91 A101 A12 A22 A32 A42 A52 ...
A62 A72 A82 A92 A102 A13 A23 A33 A43 A53 A63 A73 A83 A93 A103 A14 ...
A24 A34 A44 A54 A64 A74 A84 A94 A104 A15 A25 A35 A45 A55 A65 A75 ...
A85 A95 A105 A16 A26 A36 A46 A56 A66 A76 A86 A96 A106 A17 A27 A37 ...
A47 A57 A67 A77 A87 A97 A107 E Ta rhoA1 rhoAj F1 F2 F3 F4 F5 F6 ...
F7 F8 a Tp D1 D2 D3 D4 D5 D6 D7 D8 G1 G2 G3 G4 G5 G6 ...
u0 s0 M Tc pc rhoC T0 tauC
A = [ A11 A12 A13 A14 A15 A16 A17;
A21 A22 A23 A24 A25 A26 A27;
A31 A32 A33 A34 A35 A36 A37;
A41 A42 A43 A44 A45 A46 A47;
A51 A52 A53 A54 A55 A56 A57;
A61 A62 A63 A64 A65 A66 A67;
A71 A72 A73 A74 A75 A76 A77;
A81 A82 A83 A84 A85 A86 A87;
A91 A92 A93 A94 A95 A96 A97;
A101 A102 A103 A104 A104 A106 A107; ];
tauAj = [ 1.544912 2.5 2.5 2.5 2.5 2.5 2.5 ];
P = rho*R*T;
for j = 1:7
[C(j)] = C_fun(T,j);
[H(j)] = H_fun(rho,j);
P = P + C(j)*H(j);
end
end
If I put the right values of rho and T and start the programme, it doesn't return the right pressure (as I see in the water tables)
I think I've made a mistake in the index i,j, when I transform the C++ in code into Matlab code.
Is there someone who can help me? Thank you

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by