Function handle with integrals of multiple equations?
조회 수: 4 (최근 30일)
이전 댓글 표시
Thank you for answering
댓글 수: 6
Torsten
2020년 4월 4일
You don't need Cp since gamma=Cp/Cp-R =1-R :-)
I think you meant gamma = Cp/(Cp-R).
채택된 답변
Star Strider
2020년 4월 4일
편집: Star Strider
2020년 4월 4일
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
Star Strider
2020년 4월 6일
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
2020년 4월 7일
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 ?
추가 답변 (1개)
Ameer Hamza
2020년 4월 4일
편집: Ameer Hamza
2020년 4월 5일
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
Ameer Hamza
2020년 4월 5일
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.
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!