求转子曲线所包围的封闭区域的面积

问题

碰到这样的问题，感觉很神奇。
定子方程，短幅内摆线方程：

$$\left\{\quad \begin{array}{rl} x_1&=(R-r) \sin \tau+e \sin \left(z_2 \tau\right)-r_e \sin \theta \ y_1&=(R-r) \cos\tau-e \cos \left(z_2 \tau\right)+r_e \cos \theta \\end{array} \right.$$
与定子曲线方程共轭的转子曲线方程：
$$\left\{\quad \begin{array}{c} x_2=x_1 \cos (\varphi -\psi )-y_1 \sin (\varphi -\psi )-e \sin (\psi ) \ y_2=x_1 \sin (\varphi -\psi )+y_1 \cos (\varphi -\psi )-e \cos (\psi ) \\end{array} \right.$$

其中：
1. $R=48.78$ 为导圆半径，
2. $r=8.13$ 为滚圆半径，
3. $z_2=z_1-1$ 为转子头数。
4. $e=7.05$ 为偏心距，
5. $\theta =\tan^{-1}\dfrac{\sin\left(z_1\tau\right)}{f+\cos\left(z_1\tau\right)}-\tau$,
6. $r_e=12.6$为等距半径，
7. $\varphi=\sin^{-1}[f\sin(\theta+\tau)]$,
8. $\dfrac{\psi}{\varphi}=\dfrac{z_1}{z_1-1}$

求如下转子曲线外轮廓所包围的面积。

原理

这种情况，曲线太复杂了，符号计算的可能性太小，投入和收获之间不成比例。所以，近似的数值计算是合理的。

涉及两方面的主要问题：
1. 如何找出外轮廓；
2. 如何计算外轮廓包围的面积。

第一个问题，因为近似计算，用曲线 $\langle x(\tau),y(\tau)\rangle$ 当自己相交时对应的 $(x,y)$ 坐标相等来近似计算的。即 $\langle x(\tau_1)-x(\tau_2),y(\tau_1)-y(\tau_2)\rangle=\langle0,0\rangle$, 这是一个二元非线性方程组求根的问题，可以用牛顿法来解决。 用图示或其它方法确定出相交时曲线两部分的$\tau_1$ 和 $\tau_2$ 的大致范围可以帮助确定初值。

第二个问题，前面一篇博客刚刚提到，把要计算部分的轮廓离散化成多边形，然后用Shoelace公式计算多边形的面积即可。对应于代码中的 PolygonSignedArea 函数。

答案

直接上代码：

ClearAll["Global`*"];
R = 48.78;
r = 8.13;
z1 = R/r;
z2 = 1 - z1;
e = 7.05;
f = r/e;
re = 12.6;
θ = ArcTan[Sin[z1 τ]/(f + Cos[z1 τ])] - τ;
φ = ArcSin[f Sin[θ + τ]] - θ;
ψ = z1/(z1 - 1) φ;
curve01 = {(R - r) Sin[τ] + e Sin[z2 τ] -re Sin[θ], (R - r) Cos[τ] - e Cos[z2 τ] +re Cos[θ]} // FullSimplify;
curve02 = {curve01[[1]] Cos[φ - ψ] - curve01[[2]] Sin[φ - ψ] - e Sin[ψ], curve01[[1]] Sin[φ - ψ] +  curve01[[2]] Cos[φ - ψ] - e Cos[ψ]} //  Simplify;

base = ParametricPlot[Evaluate[curve02], {\[Tau], 0, 5 \[Pi]},Exclusions -> None, MaxRecursion -> 15, PlotPoints -> 1000];

r1 = FindRoot[ (curve02 /. τ -> x) == (curve02 /. τ -> y), {x, .5}, {y, 5.5}];
p1 = y /. FindRoot[ (ArcTan @@ (curve02 /. τ -> y)) == 
3 Pi/ 10 , {y, .55}];
top = x /. FindRoot[ curve02[[1]] /.  τ -> x , {x, 5}];
arc = Join[Table[ curve02 , {τ, top, y /. r1 , .0001}],    Table[ curve02 , {τ, x /. r1, p1 , .0001}]];

Show[base,Graphics[{{Red,  Line[{curve02 /.τ -> top, {0, 0}, curve02 /.τ -> p1}], Line@arc }}]];

PolygonSignedArea[pts_?MatrixQ] := Total[Det /@ Partition[pts, 2, 1, 1]]/2;
area = 10 PolygonSignedArea[Reverse@Join[{{0, 0}}, arc]];

Print[Style["Area is: ", Blue, 20], 
 Style[NumberForm[area, 10], Red, 23]];

Area is: 7936.859683