Function handle with integrals of multiple equations?
4 views (last 30 days)
Show older comments
Thank you for answering
6 Comments
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).
Accepted Answer
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
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
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 ?
More Answers (1)
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
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!