# Solving an equation for angle $\theta$

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.$$

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 and $$0 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:

• Note we are working with a diatomic gas; $\gamma = 3.5$. Why are you using $\gamma = 1.4$? – JD_PM Jun 15 at 10:31
• 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). – JD_PM Jun 15 at 10:35
• For a diatomic gas, $C_V=2.5R$ and $C_P=C_V+R=3.5R$. Then $\gamma=C_P/C_V=1.4$. – user5713492 Jun 15 at 15:38
• 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$... – user5713492 Jun 15 at 17:22
• 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$. – user5713492 Jun 15 at 17:26