Solving an equation for angle $\theta$

Question

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

$$dQ=\frac{1}{R}[(C_v+R)PdV+C_VVdP]\tag{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

Answer 1

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;
theta1*180/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);
theta2*180/pi
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)+ ...
    (s-1)*(r-1)/8*(theta1+sin(theta1)*cos(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)+ ...
    (s-1)*(r-1)/8*(theta2+sin(theta2)*cos(theta2)));
Q = Q2-Q1
W = pi*(r-1)*(s-1)/4
eta = W/Q
theta = linspace(0,2*pi,400);
plot((s+1)/2+(s-1)/2*cos(theta),(r+1)/2+(r-1)/2*sin(theta),'b-');
title(['Elliptical heat engine for s=' num2str(s) ', r=' num2str(r) ...
    ', \eta=' num2str(eta)]);
xlabel('V/V_1');
ylabel('P/P_1');
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);
plot(V1,c1./V1.^gamma,'r-');
c2 = ((r+1)/2+(r-1)/2*y2)*((s+1)/2+(s-1)/2*x2)^gamma;
V2 = linspace(1,s+1,400);
plot(V2,c2./V2.^gamma,'k-');
hold off;

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