# Some manipulations on (Kahraman, 1994)

[1] A. Kahraman, "Natural Modes of Planetary Gear Trains", Journal of Sound and Vibration, vol. 173, no. 1, pp. 125-130, 1994. https://doi.org/10.1006/jsvi.1994.1222.


```python
from sympy import *
init_printing()

def symb(x,y):
    return symbols('{0}_{1}'.format(x,y), type = float)
```

# Displacement vector:


```python
n = 3 # number of planets
N = n + 3 # number of degrees of freedom
crs = ['c', 'r', 's'] # carrier, ring, sun
pla = ['p{}'.format(idx + 1) for idx in range(n)] # planet
crs = crs + pla # put them together

coeff_list = symbols(crs)
c = coeff_list[0]
r = coeff_list[1]
s = coeff_list[2]

X = Matrix([symb('u', v) for v in coeff_list])
coeff_list[3:] = symbols(['p']*n)
p = coeff_list[3]
X.transpose() # Eq. (1a)
```




$\displaystyle \left[\begin{matrix}u_{c} & u_{r} & u_{s} & u_{p1} & u_{p2} & u_{p3}\end{matrix}\right]$



## Stiffness matrix:


where:
* $k_1$: mesh stiffness for the ring-planet gear pair
* $k_2$: mesh stiffness for the sun-planet gear pair
* $k_c$: carrier housing stiffness
* $k_r$: ring housing stiffness
* $k_s$: sun housing stiffness
* Diagonal 1, in red
* Diagonal 2, in grey
* Off-diagonal, in blue


```python
k_1, k_2, k_c, k_r, k_s = symbols('k_1 k_2 k_c k_r k_s', type = float)

# Diagonal 1:
K_d1 = zeros(3, 3)

K_d1[0, 0] = n*(k_1 + k_2) + k_c
K_d1[1, 1] = n* k_1 + k_r
K_d1[2, 2] = n* k_2 + k_s
K_d1[0, 1] = K_d1[1, 0] = -n*k_1
K_d1[0, 2] = K_d1[2, 0] = -n*k_2

# Diagonal 2:
K_d2 = eye(n)*(k_1 + k_2)

# Off diagonal:
K_od = zeros(n, n)
K_od[:, 0] = (k_1 - k_2)*ones(n, 1)
K_od[:, 1] = -k_1       *ones(n, 1)
K_od[:, 2] =        k_2 *ones(n, 1)

K = BlockMatrix([[K_d1, K_od.transpose()],
                 [K_od, K_d2]])
K = Matrix(K)

if(not K.is_symmetric()):
    print('error.')
K
```




$\displaystyle \left[\begin{matrix}3 k_{1} + 3 k_{2} + k_{c} & - 3 k_{1} & - 3 k_{2} & k_{1} - k_{2} & k_{1} - k_{2} & k_{1} - k_{2}\\- 3 k_{1} & 3 k_{1} + k_{r} & 0 & - k_{1} & - k_{1} & - k_{1}\\- 3 k_{2} & 0 & 3 k_{2} + k_{s} & k_{2} & k_{2} & k_{2}\\k_{1} - k_{2} & - k_{1} & k_{2} & k_{1} + k_{2} & 0 & 0\\k_{1} - k_{2} & - k_{1} & k_{2} & 0 & k_{1} + k_{2} & 0\\k_{1} - k_{2} & - k_{1} & k_{2} & 0 & 0 & k_{1} + k_{2}\end{matrix}\right]$



## Inertia matrix:


```python
M = diag(*[symb('m', v) for v in coeff_list])
M
```




$\displaystyle \left[\begin{matrix}m_{c} & 0 & 0 & 0 & 0 & 0\\0 & m_{r} & 0 & 0 & 0 & 0\\0 & 0 & m_{s} & 0 & 0 & 0\\0 & 0 & 0 & m_{p} & 0 & 0\\0 & 0 & 0 & 0 & m_{p} & 0\\0 & 0 & 0 & 0 & 0 & m_{p}\end{matrix}\right]$



## Remove ring degree of freedom


```python
X.row_del(1)
K.row_del(1)
K.col_del(1)
M.row_del(1)
M.col_del(1)
coeff_list.remove(r)
N = N - 1
```

## Coordinate transformation:

First from translational to torsional coordinates, them making the sun DOF to be the last one, making it easier to assemble a multi-stage gearbox.


```python
R_1 = diag(*[symb('r', v) for v in coeff_list])
R_1
```




$\displaystyle \left[\begin{matrix}r_{c} & 0 & 0 & 0 & 0\\0 & r_{s} & 0 & 0 & 0\\0 & 0 & r_{p} & 0 & 0\\0 & 0 & 0 & r_{p} & 0\\0 & 0 & 0 & 0 & r_{p}\end{matrix}\right]$



making the sun DOF to be the last one:


```python
N1 = N - 1

R_2 = zeros(N, N)
R_2[0, 0] = 1
R_2[1, N1] = 1
R_2[2:N, 1:N1] = eye(n)
R_2
```




$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1\\0 & 1 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 1 & 0\end{matrix}\right]$




```python
R = R_1*R_2
RMR = lambda m: transpose(R)*m*R
```

### Inertia matrix


```python
M = RMR(M)

if(not M.is_symmetric()):
    print('error in M matrix')
M
```




$\displaystyle \left[\begin{matrix}m_{c} r_{c}^{2} & 0 & 0 & 0 & 0\\0 & m_{p} r_{p}^{2} & 0 & 0 & 0\\0 & 0 & m_{p} r_{p}^{2} & 0 & 0\\0 & 0 & 0 & m_{p} r_{p}^{2} & 0\\0 & 0 & 0 & 0 & m_{s} r_{s}^{2}\end{matrix}\right]$



### Stiffness matrix


```python
K = RMR(K)

if(not K.is_symmetric()):
    print('error in K matrix')
```

The housing stiffness for both carrier and sunare null:


```python
K = K.subs([(k_c, 0), (k_s, 0)])
K
```




$\displaystyle \left[\begin{matrix}r_{c}^{2} \left(3 k_{1} + 3 k_{2}\right) & r_{c} r_{p} \left(k_{1} - k_{2}\right) & r_{c} r_{p} \left(k_{1} - k_{2}\right) & r_{c} r_{p} \left(k_{1} - k_{2}\right) & - 3 k_{2} r_{c} r_{s}\\r_{c} r_{p} \left(k_{1} - k_{2}\right) & r_{p}^{2} \left(k_{1} + k_{2}\right) & 0 & 0 & k_{2} r_{p} r_{s}\\r_{c} r_{p} \left(k_{1} - k_{2}\right) & 0 & r_{p}^{2} \left(k_{1} + k_{2}\right) & 0 & k_{2} r_{p} r_{s}\\r_{c} r_{p} \left(k_{1} - k_{2}\right) & 0 & 0 & r_{p}^{2} \left(k_{1} + k_{2}\right) & k_{2} r_{p} r_{s}\\- 3 k_{2} r_{c} r_{s} & k_{2} r_{p} r_{s} & k_{2} r_{p} r_{s} & k_{2} r_{p} r_{s} & 3 k_{2} r_{s}^{2}\end{matrix}\right]$



From that, one can write the matrices for a planetary system with $n$-planets using the following code:


```python
m_c, m_s, m_p, r_c, r_s, r_p = symbols('m_c m_s m_p r_c r_s r_p', type = float)

M_p = zeros(N, N)
M_p[0, 0] = m_c*r_c**2
M_p[N1, N1] = m_s*r_s**2
M_p[1:N1, 1:N1] = m_p*r_p**2 * eye(n)

K_p = zeros(N, N)
K_p[0, 0] = n*(k_1 + k_2)*r_c**2
K_p[N1, 0] = -n*k_2*r_s*r_c
K_p[0, N1] = -n*k_2*r_s*r_c
K_p[N1, N1] = n*k_2*r_s**2
K_p[0, 1:N1] = (k_1 - k_2)*r_c*r_p*ones(1, n)
K_p[1:N1, 0] = (k_1 - k_2)*r_c*r_p*ones(n, 1)
K_p[N1, 1:N1] = k_2*r_p*r_s*ones(1, n)
K_p[1:N1, N1] = k_2*r_p*r_s*ones(n, 1)
K_p[1:N1, 1:N1] = (k_1 + k_2)*r_p**2 * eye(n)

m_diff = abs(matrix2numpy(simplify(M_p - M))).sum()
k_diff = abs(matrix2numpy(simplify(K_p - K))).sum()

if(m_diff != 0.0):
    print('Error in M matrix.')

if(k_diff != 0.0):
    print('Error in K matrix.')
```

## Combining planet DOFs:


```python
C = zeros(N, 3)
C[   0, 0] = 1
C[  N1, 2] = 1
C[1:N1, 1] = ones(n, 1)

CMC = lambda m: transpose(C)*m*C
```

### Inertia matrix


```python
M_C = CMC(M)

if(not M_C.is_symmetric()):
    print('error in M_C matrix')
M_C
```




$\displaystyle \left[\begin{matrix}m_{c} r_{c}^{2} & 0 & 0\\0 & 3 m_{p} r_{p}^{2} & 0\\0 & 0 & m_{s} r_{s}^{2}\end{matrix}\right]$



### Stiffness matrix


```python
K_C = CMC(K)

if(not K_C.is_symmetric()):
    print('error in M_C matrix')
K_C
```




$\displaystyle \left[\begin{matrix}r_{c}^{2} \left(3 k_{1} + 3 k_{2}\right) & 3 r_{c} r_{p} \left(k_{1} - k_{2}\right) & - 3 k_{2} r_{c} r_{s}\\3 r_{c} r_{p} \left(k_{1} - k_{2}\right) & 3 r_{p}^{2} \left(k_{1} + k_{2}\right) & 3 k_{2} r_{p} r_{s}\\- 3 k_{2} r_{c} r_{s} & 3 k_{2} r_{p} r_{s} & 3 k_{2} r_{s}^{2}\end{matrix}\right]$



## Adapting it to a parallel gear set

Considering only one of the sun-planets pairs, one should change the sub-indexes in the following way:
* [p]lanet => [w]heel
* [s]un    => [p]inion;
It also necessary to remove the mesh stiffness of the ring-planet pair
### Inertia matrix


```python
k, w, p = symbols('k w p', type = float)
m_w, m_p, r_w, r_p = symbols('m_w m_p r_w r_p', type = float)

N2 = N - 2
M_par = M[N2:, N2:]
M_par = M_par.subs([(m_p, m_w), (m_s, m_p), (r_p, r_w), (r_s, r_p)]) # 
M_par
```




$\displaystyle \left[\begin{matrix}m_{w} r_{w}^{2} & 0\\0 & m_{p} r_{p}^{2}\end{matrix}\right]$



### Stiffness matrix


```python
K_par = K[N2:, N2:]

K_par = K_par.subs(k_1, 0) # ring-planet mesh stiffness
K_par = K_par.subs(k_s, 0) # sun's bearing stiffness
K_par = K_par.subs(n*k_2, k_2) # only one pair, not n
K_par = K_par.subs(k_2, k) # mesh-stiffness of the pair
K_par = K_par.subs([(r_p, r_w), (r_s, r_p)])
K_par
```




$\displaystyle \left[\begin{matrix}k r_{w}^{2} & k r_{p} r_{w}\\k r_{p} r_{w} & k r_{p}^{2}\end{matrix}\right]$



From that, one can write the matrices for a parallel system using the following code:


```python
M_p = diag(m_w*r_w**2, m_p*r_p**2)

mat_diff = abs(matrix2numpy(simplify(M_p - M_par))).sum()
if(mat_diff != 0.0):
    print('Error in M_p matrix.')

K_p       = diag(r_w**2, r_p**2)
K_p[0, 1] = r_p*r_w
K_p[1, 0] = r_p*r_w
K_p       = k*K_p

mat_diff = abs(matrix2numpy(simplify(K_p - K_par))).sum()
if(mat_diff != 0.0):
    print('Error in K_p matrix.')
```

66e8ae45f6a5669a8a16c1f9e2d5b56ff1bf17b1
ipynb
Jupyter Notebook
notes/.ipynb_checkpoints/Kahraman_1994-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
2020-10-17T13:43:01.000Z
2020-10-17T13:43:01.000Z
notes/.ipynb_checkpoints/Kahraman_1994-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
notes/.ipynb_checkpoints/Kahraman_1994-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
Qwen/Qwen-72B
1. YES
2. YES