The shape of the ellipse can be represented parametrically in terms of the angle $\theta$ around the cycle, assuming $\theta$ is measured clockwise from the point $(V_0-(V_0-V_\mathrm{max}),P_0)$:

$$P=P_0+(P_\mathrm{max}-P_0)\sin{\theta}\tag{1a}$$ $$V=V_0-(V_\mathrm{max}-V_0)\cos{\theta}\tag{1b}$$

There are two angles $\theta$ at which adiabats are tangent to the ellipse. These two angles can be obtained by substituting Eqns. 1 into Eqn. 3


with $dQ = 0$.

So what I did was deriving $P$ and $V$ wrt the angle:

$$dP = (P_{max} -P_0)\cos \theta$$

$$dV = -(V_{max} -V_0)\sin \theta$$

And plugging it into 3:

$$0=[-(C_v+R)P(V_{max} -V_0)\sin \theta+C_VV(P_{max} -P_0)\cos \theta]$$

This above equation is satisfied by two $\theta$ angles. But how to solve it?

This is what I have been told:

There are two angles $\theta$ at which adiabats are tangent to the ellipse. These two angles can be obtained by substituting Eqns. 1 into Eqn. 3 with $dQ = 0.$

More details: https://chemistry.stackexchange.com/questions/116757/how-to-get-the-efficiency-of-a-heat-engine-which-undergoes-an-elliptical-cycle/116770#116770


We assume the extremes in pressure are given by $P_1$ and $P_2=rP_1$ where $r$ is the pressure compression ratio so $$P_0=\frac{P_2+P_1}2=\frac{r+1}2P_1$$ and then $$P_{max}-P_0=P_2-P_0=rP_1-\frac{r+1}2P_1=\frac{r-1}2P_1$$ So on the ellipse $$P=\left(\frac{r+1}2+\frac{r-1}2\sin\theta\right)P_1$$ Similarly the extremes in volume are $V_1$ and $V_2=sV_1$ where $s$ is the volume compression ratio, so on the ellipse $$V=\left(\frac{s+1}2+\frac{s-1}2\cos\theta\right)V_1$$ The path is oriented such that $\theta$ decreases as the engine produces power. $P$ and $V$ then satisfy the equation $$\frac{\left(V-\frac{s+1}2V_1\right)^2}{\left(\frac{s-1}2V_1\right)^2}+\frac{\left(P-\frac{r+1}2P_1\right)^2}{\left(\frac{r-1}2P_1\right)^2}=1$$ So $$\frac{dP}{dV}=-\frac{\left(\frac{r-1}2P_1\right)^2}{\left(\frac{s-1}2V_1\right)^2}\frac{\left(V-\frac{s+1}2V_1\right)}{\left(P-\frac{r+1}2P_1\right)}$$ During adiabatic compression and expansion, $P$ and $V$ satisfy $$PV^{\gamma}=\text{constant}$$ Where $\gamma$ is the ratio of specific heats: $C_P=C_V+R=\gamma C_V$. Then $$\frac{dP}{dV}=-\gamma\frac PV=-\gamma\frac{\left(P-\frac{r+1}2P_1\right)+\frac{r+1}2P_1}{\left(V-\frac{s+1}2V_1\right)+\frac{s+1}2V_1}$$ If we equate our $2$ expressions for $\frac{dP}{dV}$ we get $$\gamma\frac{\left(P-\frac{r+1}2P_1\right)^2+\left(\frac{r+1}2P_1\right)\left(P-\frac{r+1}2P_1\right)}{\left(\frac{r-1}2P_1\right)^2}=\frac{\left(V-\frac{s+1}2V_1\right)^2+\left(\frac{s+1}2V_1\right)\left(V-\frac{s+1}2V_1\right)}{\left(\frac{s-1}2V_1\right)^2}$$ We can substitute $$\frac{V-\frac{s+1}2V_1}{\frac{s-1}2V_1}=\pm\sqrt{1-\left(\frac{P-\frac{r+1}2P_1}{\frac{r-1}2P_1}\right)^1}$$ And also let $$y=\frac{P-\frac{r+1}2P_1}{\frac{r-1}2P_1}$$ And we get $$\gamma y^2+\gamma\frac{r+1}{r-1}y=1-y^2\pm\frac{s+1}{s-1}\sqrt{1-y^2}$$ Since the adiabats are decreasing functions of $V$ we want solutions in the first and third quadrants so we take the plus sign when $y>0$ and the minus sign whe $y<0$. Squaring we get the horrible quartic equation $$\left((\gamma+1)y^2+\gamma\frac{r+1}{r-1}y-1\right)^2=\left(\frac{s+1}{s-1}\right)^2\left(1-y^2\right)$$ We can solve to get solutions $-1<y_1<0$ and $0<y_2<1$ and then find $x_1=-\sqrt{1-y_1^2}$ and $x_2=\sqrt{1-y_2^2}$. Then $\theta_1=\text{atan2}\left(y_1,x_1\right)+2\pi$, $\theta_2=\text{atan2}\left(y_2,x_2\right)$. Then we have $$dQ=\frac1R(C_V+R)P\,dV+\frac1RC_VV\,dP=\frac1{\gamma-1}\left(\gamma P\,dV+V\,dP\right)$$ $$\begin{align}Q_h&=\frac1{\gamma-1}\int_{\theta_1}^{\theta_2}\left[-\gamma\left(\frac{r+1}2P_1+\frac{r-1}2P_1\sin\theta\right)\frac{s-1}2V_1\sin\theta\right.\\ &+\left.\left(\frac{s+1}2V_1+\frac{s-1}2V_1\cos\theta\right)\frac{r-1}2P_1\cos\theta\right]d\theta\\ &=\frac{P_1V_1}{\gamma-1}\left[\gamma\left(\frac{r+1}2\right)\left(\frac{s-1}2\right)\sin\theta+\gamma\left(\frac{r-1}2\right)\left(\frac{s-1}2\right)\frac12\left(\theta+\sin\theta\cos\theta\right)\right.\\ &+\left.\left(\frac{s+1}2\right)\left(\frac{r-1}2\right)\cos\theta-\left(\frac{s-1}2\right)\left(\frac{r-1}2\right)\frac12\left(\theta-\sin\theta\cos\theta\right)\right]_{\theta_1}^{\theta_2}\end{align}$$ And $$W=\text{Area}=\pi\left(\frac{s-1}2V_1\right)\left(\frac{r-1}2P_1\right)$$ And then we can get the efficiency $$\eta=\frac W{Q_h}$$ So I worked this into a Matlab program

% ellipse.m

clear all;
close all;

r = 2; % Pressure compression ratio: P2 = r*P1
s = 3; % Volume compression ratio: V2 = s*V1
gamma = 1.4; % Ratio of specific heats
F = [(gamma+1)^2 2*gamma*(gamma+1)*(r+1)/(r-1) ...
    ((s+1)/(s-1))^2-2*(gamma+1)+gamma^2*((r+1)/(r-1))^2 ...
    -2*gamma*(r+1)/(r-1) 1-((s+1)/(s-1))^2];
yvals = roots(F);
ind1 = find(~imag(yvals) & yvals < 0 & yvals > -1 & ...
    (gamma+1)*yvals.^2+gamma*(r+1)/(r-1)*yvals-1 < 0);
y1 = yvals(ind1);
x1 = -sqrt(1-y1^2);
theta1 = atan2(y1,x1)+2*pi;
ind2 = find(~imag(yvals) & yvals > 0 & yvals < 1 & ...
    (gamma+1)*yvals.^2+gamma*(r+1)/(r-1)*yvals-1 > 0);
y2 = yvals(ind2);
x2 = sqrt(1-y2^2);
theta2 = atan2(y2,x2);
Q1 = 1/(gamma-1)*(gamma*(r+1)*(s-1)/4*cos(theta1)- ...
    gamma*(r-1)*(s-1)/8*(theta1-sin(theta1)*cos(theta1))+ ...
    (s+1)*(r-1)/4*sin(theta1)+ ...
Q2 = 1/(gamma-1)*(gamma*(r+1)*(s-1)/4*cos(theta2)- ...
    gamma*(r-1)*(s-1)/8*(theta2-sin(theta2)*cos(theta2))+ ...
    (s+1)*(r-1)/4*sin(theta2)+ ...
Q = Q2-Q1
W = pi*(r-1)*(s-1)/4
eta = W/Q
theta = linspace(0,2*pi,400);
title(['Elliptical heat engine for s=' num2str(s) ', r=' num2str(r) ...
    ', \eta=' num2str(eta)]);
hold on;
axis([0 s+1 0 r+1]);
c1 = ((r+1)/2+(r-1)/2*y1)*((s+1)/2+(s-1)/2*x1)^gamma;
V1 = linspace(1/2,s+1,400);
c2 = ((r+1)/2+(r-1)/2*y2)*((s+1)/2+(s-1)/2*x2)^gamma;
V2 = linspace(1,s+1,400);
hold off;

And plotted the ellipse and the critical adiabats in the $PV$-plane: Figure 1

  • $\begingroup$ Note we are working with a diatomic gas; $\gamma = 3.5$. Why are you using $\gamma = 1.4$? $\endgroup$ – JD_PM Jun 15 at 10:31
  • $\begingroup$ Your answer is great. However, could you add how you got $\theta_1=195.0948° $and $\theta_2=30.2732°$? I am not acquainted with solving that kind of quadratic equations (you gave two lines of explanation but I still don't know how to get the angles). $\endgroup$ – JD_PM Jun 15 at 10:35
  • 1
    $\begingroup$ For a diatomic gas, $C_V=2.5R$ and $C_P=C_V+R=3.5R$. Then $\gamma=C_P/C_V=1.4$. $\endgroup$ – user5713492 Jun 15 at 15:38
  • 1
    $\begingroup$ As for the quartic equation, I expanded mine out to the form $$ay^4+by^3+cy^2+dy+e=0$$. There is a [tag for quartic equations](math.stackexchange.com/questions/tagged/quartic-equations), you know. Since their solution implies solving a resolvent cubic, you probably also need the [tag for cubic equations](math.stackexchange.com/questions/tagged/cubic-equations). Rather than do all that I just used Matlab's [roots](mathworks.com/help/matlab/ref/roots.html) function. Then knowing $y=\sin\theta$ and $x=\pm\sqrt{1-y^2}=\cos\theta$ has the same sign as $y$... $\endgroup$ – user5713492 Jun 15 at 17:22
  • $\begingroup$ Then Matlab has a function that is useful for rectangular to polar conversions, atan2 which got me $\theta$. Since $\theta$ decreases as the engine produces work I needed to add $2\pi$ to my value for $\theta_1$ to make it bigger than $\theta_2$. $\endgroup$ – user5713492 Jun 15 at 17:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.