参考文献：Bell I H. Theoretical and experimental analysis of liquid flooded compression in scroll compressors[J]. Dissertations & Theses - Gradworks, 2011.
- 图中φ1和φ2的位置错了，应该是逆时针布置才对，遵循右手定律
- 展开角逆时针为正
- 所以从表达式中我们可以看到动涡盘的第一项都是负的，至于为什么第二项还有一个Δθ呢？这是因为保证在一开始装配的时候就能保持啮合的状态
- 所以应该则么求Δθ呢？可以在初始状态的时候，使θ=0，动静涡盘在φie下啮合
- 求解关于初始装配条件下待定系数后的动涡盘等效曲柄转角 theta_m = geo.phi_fie - theta + 3.0*pi/2.0 = geo.phi_fie - theta - pi/2.0
- 其实本质上来讲排气角的条件是啮合角不得小于展开开始角，用这个原则来判断，静涡盘内壁面啮合角φfic=φe-θ-2πNc;动涡盘外壁面啮合角φfic=φe-π-θ-2πNc;
- 如果采用静涡盘内壁面啮合角φfic = φe-θ-2πNc，则由于φfic > φfis，即φfis < φe-θ-2πNc
- 如果采用动涡盘内壁面啮合角φoic = φe-π-θ-2πNc，则由于φoic > φois，即φois < φe-π-θ-2πNc
- floor()函数代表向下取整
- 第一种定义：固定在动涡盘外壁面，当地展开角固定为为φie-π
- 第二种定义：固定在动涡盘外壁面，该点在静涡盘展开结束角的法线与动涡盘外壁面的交点处
- 第三种定义：固定在动涡盘外壁面，该点在静涡盘内壁面展开结束角的终点跟静涡盘基圆圆心的连线，与动涡盘外壁面的交点处
- 公转半径只是涡盘基圆半径和内外表面发生角的函数
- t The thickness of the scroll 涡旋的厚度 t=ts = rb*(phi_i0 - phi_o0)，涡旋盘厚度通常为5毫米量级。
- 密度偏移，要注意的一点是求质心坐标的时候，分子里面的坐标要求是微元的重心坐标，对于矩形微元而言重心坐标即在其中心，但是在扇形微元的时候，重心在其半径的2/3处的地方，造成了所谓的“密度偏移”
- 吸入腔内部涡旋面积体积质心计算
- 吸入腔总的面积体积质心计算
- 压缩机排量定义为动涡盘转一圈后的瞬间吸入腔的总体积，但是这个我没有推过，不知道排量的公式对还是错
- 中间是叉乘，rO is a vector from (xO, yO)，φie-pi/2之前也讲过，是动涡盘和曲柄销的安装角度
- 吸气腔水平面力和力矩的计算 但是在计算力矩的时候有点奇怪，按道理来讲，应该有一个因子是关于回转半径ro的，但是这里没有，只有rb平方。
- 其实就是在花篮中涡旋盘以外的区域
- 吸气面积腔，其实这个腔并不是很重要，因为以来这里是用来引导进气，二来这个区域一般有有其他用于加工过程的金属占据着空间
- 可以大概认为可以使用对称面积*2的方法，但是这个方法不是绝对精确的，因为吸入破入角对于静涡盘和动涡盘来讲是不相同的
- 文献中的方法，我试过是正确的
- 我的方法，形式跟文献的基本一致
- 关于压缩腔的质心的证明实在太麻烦了，这里我没有去验证
- 这个没什么好说的，跟吸入腔的计算方法类似
- 一般过了排气角π/2之后才使得d1和dd压力相等
- d1和dd腔的分界线是由动盘的φo和静盘的φos+π连线。目的是保证在排气角处分界线的长度为0
- 面积求解的套路一般都是把动涡盘的参考系转化为静涡盘的参考系
- Ad1 = AO - A1
- A1 = A1a + A1b + A1c + A1d A1将关于动涡盘的面积转化到静涡盘上
- Ad1 = AO - (A1a + A1b + A1c + A1d)
- 转化为静涡盘的参考系，排气角看起来为什么要加上π，但是你从动涡盘的角度去看，转化为动涡盘的参考系展开角是需要减去π的，此时就剩下φos - φo0了
- 动涡盘上的面积通过旋转和平移搬到静涡盘上，展开角不变
- 动涡盘上的面积通过旋转和平移搬到静涡盘上，展开角不变，运用圆渐开线基圆三角形面积求解的结论可以直接得到答案
- 运用圆渐开线基圆三角形面积求解的结论可以直接得到答案
- V1d,d1这个三角形比较麻烦，因为它的一个顶点在静涡盘上，另一个顶点在动涡盘上，算起来比较麻烦，这里我就不推导了
- A1 = A1a + A1b + A1c + A1d
- Ad1 = AO - A1
- 还是以前的套路
- dd腔的定义比较麻烦，因为它的形成曲线有不同的类型
- 从φ0~φs这一段涡旋是用来形成排气区域的组合曲线的
- 组合曲线满足四个条件
- 1，曲线必须在φis处跟涡旋内壁相切；
- 2，曲线必须在φos处跟涡旋外壁相切；
- 3，曲线在运行过程中不会使得动静涡盘相碰；
- 4，曲线必须经过涡旋的内外φs点
- 曲线类型
- 第一种组合曲线类型——弧线-直线-弧线 ra2是连接外涡旋壁面的圆弧
- 第二种组合曲线类型——弧线-弧线
- 组合曲线圆心坐标
- 组合曲线圆心半径
- 如果要保持啮合状态的话必须要满足的条件：φis = φos + π 否则后面还没有到达排气角，动涡盘就会被抬起，动静涡盘就不会再贴合
- 组合曲线啮合的第一种情况：弧线1的圆心在静盘基圆圆心，弧线2的圆心在动涡盘基圆圆心，此时有ra1 = ra2 + ro
- 组合曲线啮合的第二种情况：弧线1的圆心偏离静盘基圆圆心，弧线2的圆心在动涡盘基圆圆心，此时有ra1 = ra2 + ro
- 可以求出公切线的长度和公切点的坐标
- dd腔质心
- dd腔面积
Add = 2(AOa;dd + AOb;dd + AOc;dd + AIa;dd + AIb;dd)
- VOc;dd是由于φi ≠ φos + π，一个不完美的啮合所造成的间隙
- ta1,1指的是ACR1的公切点所在的转角
- 现有泄漏长度计算的主要问题是，它们对于处理任意数量的压缩室不是足够通用的。
- 径向泄漏 = 泄漏弧长 * 径向gap宽度�radial
- 泄漏面积的定义依赖于给定曲柄转角下压缩腔的对数
- 静涡盘展开角在φie~φie-θd的范围内存在啮合点时才会有压缩腔，压缩腔的对数Nc,max=floor[(φie-π-φis)/2π]，或者=啮合点个数/2-1，这是应为要形成压缩腔至少头跟尾各需要啮合点
- 弧线长的求解
- 三种情况
- 吸气腔和压缩腔
- 压缩腔与压缩腔
- 压缩腔和排气腔
- 计算的时候需要知道动静涡盘的径向间隙，然后再乘以齿高
- 包含展开角度的简单结构
- phi: 展开角φ
- theta: 曲轴转角
- double phi_max: 沿渐开线的最大展开角
- double phi_min: 沿渐开线的最小展开角
- double phi_0： 沿着这个渐开线发生角
- double dphi_max_dtheta: 沿渐开线的最大展开角对曲轴转角θ的导数
- double dphi_min_dtheta: 沿渐开线的最小展开角对曲轴转角θ的导数
- involute_index(类名) involute: 渐开线的渐开线指数
//一种通用的涡旋角度结构体
cdef class CVInvolute:
def __init__(self):
pass
// 渐开线类别
cdef enum involute_index:
INVOLUTE_FI
INVOLUTE_FO
INVOLUTE_OI
INVOLUTE_OO
// 渐开线类别返回字符串
cpdef bytes involute_index_to_key(int index):
"""
Return the string associated with a given index from the common_scroll_geo.involute_index enumeration
"""
if index == INVOLUTE_FI:
return bytes('fi')
elif index == INVOLUTE_FO:
return bytes('fo')
elif index == INVOLUTE_OI:
return bytes('oi')
elif index == INVOLUTE_OO:
return bytes('oo')
else:
return bytes('')
- CVInvolute Inner: 该腔室的内部渐开线
- CVInvolute Outer: 该腔室的外部渐开线
- Boolean has_line_1: 是否存在line #1
- Boolean has_line_2: 是否存在line #2
//申明内部渐开线和外部渐开线两种角度结构
cdef class CVInvolutes:
def __init__(self):
self.Inner = CVInvolute.__new__(CVInvolute)
self.Outer = CVInvolute.__new__(CVInvolute)
//打印用
def __repr__(self):
s = ''
s += "Outer.involute = {i:s}\n".format(i=involute_index_to_key(self.Outer.involute))
s += "Outer.phi_0 = {i:g}\n".format(i=self.Outer.phi_0)
s += "Outer.phi_max = {i:g}\n".format(i=self.Outer.phi_max)
s += "Outer.phi_min = {i:g}\n".format(i=self.Outer.phi_min)
s += "Inner.involute = {i:s}\n".format(i=involute_index_to_key(self.Inner.involute))
s += "Inner.phi_0 = {i:g}\n".format(i=self.Inner.phi_0)
s += "Inner.phi_max = {i:g}\n".format(i=self.Inner.phi_max)
s += "Inner.phi_min = {i:g}\n".format(i=self.Inner.phi_min)
s += "has_line_1 = {i:g}\n".format(i=self.has_line_1)
s += "has_line_2 = {i:g}".format(i=self.has_line_2)
return s
- h 齿高
- ro orbiting radius 公转半径 公转半径只是涡盘基圆半径和内外表面发生角的函数
- rb the radius of the base circle 基圆半径
- t The thickness of the scroll 涡旋的厚度 t=ts = rb*(phi_i0 - phi_o0)，涡旋盘厚度通常为5毫米量级。
- phi_fie Inner Ending Angle 静涡盘内部渐开线结束角
- phi_fis Inner Starting Angle 静涡盘内部渐开线开始角
- phi_fi0 Inner Initial Angle 静涡盘内部渐开线发生角
- phi_foe Inner Ending Angle 静涡盘外部渐开线结束角
- phi_fos Inner Starting Angle 静涡盘外部渐开线开始角
- phi_fo0 Inner Initial Angle 静涡盘外部渐开线发生角
- phi_oie Inner Ending Angle 动涡盘内部渐开线结束角
- phi_ois Inner Starting Angle 动涡盘内部渐开线开始角
- phi_oi0 Inner Initial Angle 动涡盘内部渐开线发生角
- phi_ooe Inner Ending Angle 动涡盘外部渐开线结束角
- phi_oos Inner Starting Angle 动涡盘外部渐开线开始角
- phi_oo0 Inner Initial Angle 动涡盘外部渐开线发生角
- phi_ie_offset 初始值为0
- copy_inplace(self, geoVals target) 结构性复制
- is_symmetric(self) → bool 动外部渐开线所有角度都相等的时候返回true
- val_if_symmetric(self, double val) → double 如果is_symmetric则返回数值，否则返回数值错误
- xvec_disc_port 属于numpy.ndarray类型
- yvec_disc_port 属于numpy.ndarray类型
- rebuild_geoVals
- setattr(x, 'y', v) is equivalent to ``x.y = v''
- getattr(x, 'y') is equivalent to x.y
- repr 打印所有属性
#This is a list of all the members in geoVals
geoValsvarlist=['h','ro','rb','t',
'phi_fi0','phi_fis','phi_fie',
'phi_fo0','phi_fos','phi_foe',
'phi_oi0','phi_ois','phi_oie',
'phi_oo0','phi_oos','phi_ooe',
'xa_arc1','ya_arc1','ra_arc1','t1_arc1','t2_arc1',
'xa_arc2','ya_arc2','ra_arc2','t1_arc2','t2_arc2',
'b_line', 't1_line', 't2_line', 'm_line',
'x0_wall','y0_wall','r_wall',
'delta_radial', 'delta_flank',
'phi_ie_offset','delta_suction_offset',
'cx_scroll','cy_scroll','V_scroll','Vremove']
//外部渐开线所有角度都相等的时候返回true
cpdef bint is_symmetric(self):
"""
Returns true if all the angles for the fixed scroll are the same as for the orbiting scroll
"""
return (abs(self.phi_fi0-self.phi_oi0) < 1e-14
and abs(self.phi_fis-self.phi_ois) < 1e-14
and abs(self.phi_fie-self.phi_oie) < 1e-14
and abs(self.phi_fo0-self.phi_oo0) < 1e-14
and abs(self.phi_fos-self.phi_oos) < 1e-14
and abs(self.phi_foe-self.phi_ooe) < 1e-14
)
- 面积的不定积分项
cpdef double Gr(double phi, geoVals geo, double theta, int inv):
"""
The antiderivative of the area integration term, where
.. math::
Gr \\equiv \\int\\left[\\left(-y\\frac{dx(\\phi)}{d\\phi}+x\\frac{dy(\\phi)}{d\\phi}\\right)d\\phi\\right]
"""
//求解关于初始装配条件下待定系数后的动涡盘等效曲柄转角
theta_m = geo.phi_fie - theta + 3.0*pi/2.0
if inv == INVOLUTE_FI:
return phi*geo.rb**2*(phi**2 - 3*phi*geo.phi_fi0 + 3*geo.phi_fi0**2)/3
elif inv == INVOLUTE_FO:
return phi*geo.rb**2*(phi**2 - 3*phi*geo.phi_fo0 + 3*geo.phi_fo0**2)/3
elif inv == INVOLUTE_OI:
return geo.rb*(phi**3*geo.rb - 3*phi**2*geo.phi_oi0*geo.rb
+ 3*phi*geo.phi_oi0**2*geo.rb
+ 3*(phi-geo.phi_oi0)*geo.ro*cos(phi - theta_m)
- 3*geo.ro*sin(phi - theta_m))/3
elif inv == INVOLUTE_OO:
return geo.rb*(phi**3*geo.rb - 3*phi**2*geo.phi_oo0*geo.rb
+ 3*phi*geo.phi_oo0**2*geo.rb
+ 3*(phi-geo.phi_oo0)*geo.ro*cos(phi - theta_m)
- 3*geo.ro*sin(phi - theta_m))/3
- 面积积分的对展开角φ的偏导数
cpdef double dGr_dphi(double phi, geoVals geo, double theta, int inv):
"""
The partial derivative of Gr with respect to phi with theta held constant
"""
THETA = geo.phi_fie - theta - pi/2.0
if inv == INVOLUTE_FI:
return geo.rb**2*(phi - geo.phi_fi0)**2
elif inv == INVOLUTE_FO:
return geo.rb**2*(phi - geo.phi_fo0)**2
elif inv == INVOLUTE_OI:
return geo.rb*(geo.rb*(phi - geo.phi_oi0)**2 + (phi- geo.phi_oi0)*geo.ro*sin(THETA - phi))
elif inv == INVOLUTE_OO:
return geo.rb*(geo.rb*(phi - geo.phi_oo0)**2 + (phi- geo.phi_oo0)*geo.ro*sin(THETA - phi))
- 面积积分的对曲轴转角θ的偏导数
cpdef double dGr_dtheta(double phi, geoVals geo, double theta, int inv):
"""
The partial derivative of Gr with respect to theta with phi held constant
"""
THETA = geo.phi_fie - theta - pi/2.0
if inv == INVOLUTE_FI or inv == INVOLUTE_FO:
return 0.0
elif inv == INVOLUTE_OI:
return geo.rb*geo.ro*((phi - geo.phi_oi0)*sin(THETA - phi) - cos(THETA - phi))
elif inv == INVOLUTE_OO:
return geo.rb*geo.ro*((phi - geo.phi_oo0)*sin(THETA - phi) - cos(THETA - phi))
- 与涡旋压缩机腔室传热计算所需计算有关的角度的结构
- double phi_1_i 内部渐开线外壁最大展开角
- double phi_2_i 内部渐开线外壁最小展开角
- double phi_1_o 外部渐开线内壁最大展开角
- double phi_2_o 外部渐开线内壁最小展开角
- double phi_i0 内部渐开线外壁原始展开角
- double phi_o0 外部渐开线内壁原始展开角
cdef class HTAnglesClass:
def __repr__(self):
s = ''
for k in ['phi_1_i','phi_2_i','phi_1_o','phi_2_o','phi_i0','phi_o0']:
s += k + ' : ' + str(getattr(self,k)) + '\n'
return s
- 以一种通用的方式评估V和dV/dtheta[m^3/radian]
- cpdef VdVstruct VdV(double theta, geoVals geo, CVInvolutes inv)
cpdef VdVstruct VdV(double theta, geoVals geo, CVInvolutes inv):
"""
Evaluate V and dV/dtheta in a generalized manner for a chamber
"""
cdef double A_i, A_o, A_line_1 = 0, A_line_2 = 0, x_1, x_2, y_1, y_2, dx_1_dphi=0, dy_1_dphi=0, dx_2_dphi=0, dy_2_dphi=0
cdef double dA_line_1_dtheta = 0, dA_line_2_dtheta = 0, dx_1_dtheta, dy_1_dtheta, dx_2_dtheta, dy_2_dtheta
## ------------------------ VOLUME -------------------------------
A_i = 0.5*(Gr(inv.Outer.phi_max, geo, theta, inv.Outer.involute) - Gr(inv.Outer.phi_min, geo, theta, inv.Outer.involute))
if inv.has_line_1:
_coords_inv_d_int(inv.Outer.phi_max, geo, theta, inv.Outer.involute, &x_1, &y_1)
_coords_inv_d_int(inv.Inner.phi_max, geo, theta, inv.Inner.involute, &x_2, &y_2)
A_line_1 = 0.5*(x_1*y_2 - x_2*y_1)
A_o = 0.5*(Gr(inv.Inner.phi_min, geo, theta, inv.Inner.involute) - Gr(inv.Inner.phi_max, geo, theta, inv.Inner.involute))
if inv.has_line_2:
_coords_inv_d_int(inv.Inner.phi_min, geo, theta, inv.Inner.involute, &x_1, &y_1)
_coords_inv_d_int(inv.Outer.phi_min, geo, theta, inv.Outer.involute, &x_2, &y_2)
A_line_2 = 0.5*(x_1*y_2 - x_2*y_1)
V = geo.h*(A_i + A_line_1 + A_o + A_line_2)
## ------------------------ DERIVATIVE -------------------------------
dA_i_dtheta = 0.5*(dGr_dphi(inv.Outer.phi_max, geo, theta, inv.Outer.involute)*inv.Outer.dphi_max_dtheta
+dGr_dtheta(inv.Outer.phi_max, geo, theta, inv.Outer.involute)
-dGr_dphi(inv.Outer.phi_min, geo, theta, inv.Outer.involute)*inv.Outer.dphi_min_dtheta
-dGr_dtheta(inv.Outer.phi_min, geo, theta, inv.Outer.involute)
)
if inv.has_line_1:
coords_inv_dtheta(inv.Outer.phi_max, geo, theta, inv.Outer.involute, &dx_1_dtheta, &dy_1_dtheta)
coords_inv_dtheta(inv.Inner.phi_max, geo, theta, inv.Inner.involute, &dx_2_dtheta, &dy_2_dtheta)
_dcoords_inv_dphi_int(inv.Outer.phi_max, geo, theta, inv.Outer.involute, &dx_1_dphi, &dy_1_dphi)
dx_1_dtheta += dx_1_dphi*inv.Outer.dphi_max_dtheta
dy_1_dtheta += dy_1_dphi*inv.Outer.dphi_max_dtheta
_dcoords_inv_dphi_int(inv.Inner.phi_max, geo, theta, inv.Inner.involute, &dx_2_dphi, &dy_2_dphi)
dx_2_dtheta += dx_2_dphi*inv.Inner.dphi_max_dtheta
dy_2_dtheta += dy_2_dphi*inv.Inner.dphi_max_dtheta
dA_line_1_dtheta = 0.5*(x_1*dy_2_dtheta + y_2*dx_1_dtheta - x_2*dy_1_dtheta - y_1*dx_2_dtheta)
dA_o_dtheta = 0.5*(dGr_dphi(inv.Inner.phi_min, geo, theta, inv.Inner.involute)*inv.Inner.dphi_min_dtheta
+dGr_dtheta(inv.Inner.phi_min, geo, theta, inv.Inner.involute)
-dGr_dphi(inv.Inner.phi_max, geo, theta, inv.Inner.involute)*inv.Inner.dphi_max_dtheta
-dGr_dtheta(inv.Inner.phi_max, geo, theta, inv.Inner.involute)
)
if inv.has_line_2:
coords_inv_dtheta(inv.Inner.phi_min, geo, theta, inv.Inner.involute, &dx_1_dtheta, &dy_1_dtheta)
coords_inv_dtheta(inv.Outer.phi_min, geo, theta, inv.Outer.involute, &dx_2_dtheta, &dy_2_dtheta)
_dcoords_inv_dphi_int(inv.Inner.phi_min, geo, theta, inv.Inner.involute, &dx_1_dphi, &dy_1_dphi)
dx_1_dtheta += dx_1_dphi*inv.Inner.dphi_min_dtheta
dy_1_dtheta += dy_1_dphi*inv.Inner.dphi_min_dtheta
_dcoords_inv_dphi_int(inv.Outer.phi_min, geo, theta, inv.Outer.involute, &dx_2_dphi, &dy_2_dphi)
dx_2_dtheta += dx_2_dphi*inv.Outer.dphi_min_dtheta
dy_2_dtheta += dy_2_dphi*inv.Outer.dphi_min_dtheta
dA_line_2_dtheta = 0.5*(x_1*dy_2_dtheta + y_2*dx_1_dtheta - x_2*dy_1_dtheta - y_1*dx_2_dtheta)
dV = geo.h*(dA_i_dtheta + dA_line_1_dtheta + dA_o_dtheta + dA_line_2_dtheta)
cdef VdVstruct VdV = VdVstruct.__new__(VdVstruct)
VdV.V = V
VdV.dV = dV
return VdV
- 对应于沿着渐开线（外部渐开线盘内部[fi]，外部渐开线盘外部渐开线[fo]，内部渐开线外部渐开线[oo]和内部渐开线内部渐开线[oi]）的点的渐开线角度
- coords_inv(phi, geoVals geo, double theta, flag='fi') → tuple
- phi_vec – 渐开线角度的一维数字数组或双精度矢量
- geo – geoVals class scroll compressor geometry
- theta – 单精度曲轴转角θ 0~2π
- flag – 渐开线字符串flag： ‘fi’,’fo’,’oi’,’oo’
- return - 涡旋盘坐标元组(数组) 如：(x,y)
def coords_inv(phi,geo,theta,flag="fi")
if type(phi) is np.ndarray:
return _coords_inv_np(phi,geo,theta,flag)
else:
return _coords_inv_d(phi,geo,theta,flag)
//
cpdef tuple _coords_inv_np(np.ndarray[np.float_t] phi, geoVals geo,double theta, flag=""):
"""
Internal function that does the calculation if phi is definitely a 1D numpy vector
"""
rb = geo.rb
ro = rb*(pi - geo.phi_fi0 + geo.phi_oo0)
om = geo.phi_fie - theta + 3.0*pi/2.0
if flag=="fi":
x = rb*np.cos(phi)+rb*(phi-geo.phi_fi0)*np.sin(phi)
y = rb*np.sin(phi)-rb*(phi-geo.phi_fi0)*np.cos(phi)
elif flag=="fo":
x = rb*np.cos(phi)+rb*(phi-geo.phi_fo0)*np.sin(phi)
y = rb*np.sin(phi)-rb*(phi-geo.phi_fo0)*np.cos(phi)
elif flag=="oi":
x = -rb*np.cos(phi)-rb*(phi-geo.phi_oi0)*np.sin(phi)+ro*np.cos(om)
y = -rb*np.sin(phi)+rb*(phi-geo.phi_oi0)*np.cos(phi)+ro*np.sin(om)
elif flag=="oo":
x = -rb*np.cos(phi)-rb*(phi-geo.phi_oo0)*np.sin(phi)+ro*np.cos(om)
y = -rb*np.sin(phi)+rb*(phi-geo.phi_oo0)*np.cos(phi)+ro*np.sin(om)
else:
raise ValueError('flag not valid')
return (x,y)
//
cpdef tuple _coords_inv_d(double phi, geoVals geo,double theta, flag=""):
"""
Internal function that does the calculation if phi is a double variable
"""
rb = geo.rb
ro = rb*(pi - geo.phi_fi0 + geo.phi_oo0)
om = geo.phi_fie - theta + 3.0*pi/2.0
if flag=="fi":
x = rb*cos(phi)+rb*(phi-geo.phi_fi0)*sin(phi)
y = rb*sin(phi)-rb*(phi-geo.phi_fi0)*cos(phi)
elif flag=="fo":
x = rb*cos(phi)+rb*(phi-geo.phi_fo0)*sin(phi)
y = rb*sin(phi)-rb*(phi-geo.phi_fo0)*cos(phi)
elif flag=="oi":
x = -rb*cos(phi)-rb*(phi-geo.phi_oi0)*sin(phi)+ro*cos(om)
y = -rb*sin(phi)+rb*(phi-geo.phi_oi0)*cos(phi)+ro*sin(om)
elif flag=="oo":
x = -rb*cos(phi)-rb*(phi-geo.phi_oo0)*sin(phi)+ro*cos(om)
y = -rb*sin(phi)+rb*(phi-geo.phi_oo0)*cos(phi)+ro*sin(om)
else:
raise ValueError('flag not valid')
return (x,y)