Function handle with integrals of multiple equations?

1 view (last 30 days)
Thank you for answering
  6 Comments
Torsten
Torsten on 4 Apr 2020
You don't need Cp since gamma=Cp/Cp-R =1-R :-)
I think you meant gamma = Cp/(Cp-R).
Torsten
Torsten on 4 Apr 2020
gamma=Cp/(Cp-1) ?
I thought gamma=Cp/(Cp-R) ...

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 4 Apr 2020
Edited: Star Strider on 4 Apr 2020
I get different result with a strictly numeric version:
V1 = 230;
P1 = 2.7;
T1 = 300;
V2 = 30;
A = -0.703029;
B = 108.4773;
C = -42.52157;
D = 5.862788;
E = 0.678565;
R=8;
Cp = @(t) A+B*t+C*(t.^2)+D*(t.^3)+E./(t.^2);
gamma = @(t) Cp(t)/(Cp(t)-R);
T = @(t) 1000*t;
Eq1 = @(P2,T2) log(P2)-log(P1) - integral(@(t)(gamma(t)/(gamma(t)-1))./T(t), T1, T2);
Eq2 = @(T2) log(V2)-log(V1) - integral(@(t)(gamma(t)/(gamma(t)-1))./T(t), T1, T2);
B0 = randi(100,2,1);
B = fsolve(@(b) [Eq1(b(1),b(2)); Eq2(b(2))], B0); % Choose Appropriate Valuew For ‘B0’ To Get The Correct Result
fprintf(1, '\nP2 = %8.4f\nT2 = %8.4f\n', B)
producing:
P2 = 0.3522
T2 = 299.9684
EDIT — (4 Apr 2020 at 18:23)
V1 = 230;
P1 = 2.7;
T1 = 300;
V2 = 30;
A = -0.703029;
B = 108.4773;
C = -42.52157;
D = 5.862788;
E = 0.678565;
R=8;
Cp = @(t) A+B*t+C*(t.^2)+D*(t.^3)+E./(t.^2);
gamma = @(t) Cp(t)./(Cp(t)-R);
T = @(t) 1000*t;
Vv = V1:-0.5:V2;
B0 = randi(500,1,2);
for k = 1:numel(Vv)
Eq1 = @(P2,T2) log(P2)-log(P1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
Eq2 = @(T2) log(Vv(k))-log(V1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
Prms(k,:) = fsolve(@(b) [Eq1(b(1),b(2)); Eq2(b(2))], B0); % Choose Appropriate Valuew For ‘B0’ To Get The Correct Result
Prms = abs(Prms);
end
Out = table(Vv.', Prms(:,1), Prms(:,2), 'VariableNames',{'V','P2','T2'})
Sample output:
Out =
401×3 table
V P2 T2
_____ _______ ______
230 2.7 300
229.5 2.6941 300
229 2.6883 300
228.5 2.6824 300
228 2.6765 300
. . .
31 0.36391 299.97
30.5 0.35804 299.97
30 0.35217 299.97
  14 Comments
Star Strider
Star Strider on 6 Apr 2020
Should this:
P2 = P1*(T2/T2)^(gamma/(gamma-1))
be ‘T1/T2’ or ‘T2/T1’? I doubt that would work anyway, because ‘gamma’ needs to be a function of ‘T’ for the rest of it to work.
I was an undergraduate Chemistry major (long ago), so encountered the gas laws and the approximations to deal with non-ideal gases, however I admit to being lost with this latest set of equations. The problem is that ‘gamma’ is a function only of ‘T’, and there is no existing relation that would work in the first equation, making it integrable as a function of ‘P’. In the absence of such a relation, I doubt that it is currently possible to procede further with this.
Torsten
Torsten on 7 Apr 2020
Isn't it true that gamma(P) in the first equation has to be evaluated at the value of T which satisfies the second equation with P2 being replaced by P and T2 being replaced by T ?

Sign in to comment.

More Answers (1)

Ameer Hamza
Ameer Hamza on 4 Apr 2020
Edited: Ameer Hamza on 5 Apr 2020
This solves the equations using symbolic equations
syms V P T P2 T2
A= -0.703029;
B= 108.4773;
C= -42.52157;
D= 5.862788;
E= 0.678565;
P1= 2.7;
T1= 300;
V2= 30;
R=8;
t = T/1000;
Cp= A+B*t+C*(t^2)+D*(t^3)+E/(t^2);
gamma = Cp/(Cp-R);
V_val = 230:-0.5:V2;
P2_sol = zeros(size(V_val));
T2_sol = zeros(size(V_val));
gamma_vec = zeros(size(V_val));
for i=1:numel(V_val)
V1 = V_val(i);
eq1_lhs = int(1/P, P1, P2);
eq1_rhs = int(gamma/(gamma-1)/T, T1, T2);
eq2_lhs = int(1/V, V1, V2);
eq2_rhs = -int(1/(gamma-1)/T, T1, T2);
sol = vpasolve([eq1_lhs==eq1_rhs, eq2_lhs==eq2_rhs], [P2 T2]);
P2_sol(i) = sol.P2; % solution for P2
T2_sol(i) = sol.T2; % solution for T2
gamma_vec(i) = subs(gamma, sol.T2);
end
%%
T = table(V_val', real(P2_sol)', T2_sol', gamma_vec', 'VariableNames', {'V1', 'P2', 'T2', 'gamma'});
Since this is using the symbolic toolbox, so the speed of execution can be slow. It will take a few minutes to finish.
[Note] As you mentioned in comment to Star Strider's comment, the volume is decreasing, but remember we start with V1 = 230, and end at V2 = 30, so in that case, we will maximum change in volume, and hence the maximum change in temperature and pressure. Now suppose we start at V1 = 150 and end at V1 = 30, the difference in volume is small, and therefore the change in temperature and pressure will also be small. I hope this clearify the confusion about decreasing values of P2 and T2 when we decrease V1.
  10 Comments
Deema Khunda
Deema Khunda on 5 Apr 2020
Hey ,
I got an error in line 15 as
Unrecognized function or variable 'V_val'.
Error in combustion_trial (line 15)
P2_sol = zeros(size(V_val));
I think line line 15 V_val needs to be updated to V1_val as its defined earlier in the code as V1_val , I have updated it to the following :
syms V P T P2 T2
A= -0.703029;
B= 108.4773;
C= -42.52157;
D= 5.862788;
E= 0.678565;
P1= 2.7;
T1= 300;
V2= 30;
R=8;
t = T/1000;
Cp= A+B*t+C*(t^2)+D*(t^3)+E/(t^2);
gamma = Cp/(Cp-R);
V_val = 230:-0.5:V2;
P2_sol = zeros(size(V_val));
T2_sol = zeros(size(V_val));
for i=1:numel(V_val)
V1 = V_val(i);
eq1_lhs = int(1/P, P1, P2);
eq1_rhs = int(gamma/(gamma-1)/T, T1, T2);
eq2_lhs = int(1/V, V1, V2);
eq2_rhs = -int(1/(gamma-1)/T, T1, T2);
sol = vpasolve([eq1_lhs==eq1_rhs, eq2_lhs==eq2_rhs], [P2 T2]);
P2_sol(i) = sol.P2; % solution for P2
T2_sol(i) = sol.T2; % solution for T2
end
%%
T = table(V_val', P2_sol', T2_sol', 'VariableNames', {'V1', 'P2', 'T2'});
Can I ask if I could add gamma to the table to see how its changing with temperature?
Thank you
Ameer Hamza
Ameer Hamza on 5 Apr 2020
Deema, thats correct. That was a mistake in my code. Please check the updated answer. The value of gamma at T=T2 is also added to the table.

Sign in to comment.

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!