12 KiB
Author: Antoine Foucault-Castelli Date: 17/03/2026 Bibliography: @karmiFiabiliteOptimisationSystemes @ghorbelEffectBrakeLocation2020a @zhengAnalyticalApproachMesh2022 Version: 4.0
Gearbox simulation
Structure
main.m
func.m
time_sampling
inertia
gearbox_geometry
mesh_stiffness
wind_excitation
mass_matrix
stiffness_matrix
stiffness_update
damping_matrix
forces_matrix
initialization
newmark
parameters.m
display.m
disp_xy(t,x,i)
speed_xy(t,v,i)
acc_xy(t,v,i)
disp_theta(t,x,i)
speed_theta(t,v,i)
acc_theta(t,a,i)
mesh_stiffness(t,k_var,k_mean,i)
taero(t,T_aero,i)
rms_on_x(x)
spectral_on_x(x,i,tf,dt)
plot_x(t,x,i)
only_display.m
display_update.m
main.m is the launcher for the full code compilation
func is the function class, containing every function used in the code
parameters is the model parameters document
display is the displayer for figures
only_display to skip non-needed full compilation
display_update to update a single graph if needed
How does it work ?
For Matlab under windows : open matlab and launch main.m
Warning
remove
graphics_toolkit("gnuplot")which is only useful for Octave visualization
For Octave/Linux :
cd [path_to_cloned]/code/versions/[current_version]/
octave
>>> run main.m
>>>
Informations
- Verified code with base code and paper
- Implemented wind excitation
- Implemented Ghorbel parameters
- [~] Add varying output torque depending on number of tooth in contact ?
- [~] add brake ?
- [~] Variation of C ?
- ODE45 ?
- Displays from last version (in class)
- New mesh stiffness variation model (trapezoidal)
- Print in pdf
- Errors/defects implementation
- profile error
- assembly defect (change on gear stiffness)
- README
- Maths
- Images
- Source(s)
- Random models implementation
- Display
- time response on a2
- speed response
- x movement (vs y movement ?)
- RMS STD Kurtosis
- Tests
- different speed
- different blade radii
Maths
The model is represented by two 4DDL shafts in a gear connection : x and y are free at the pinion, \theta is free at both ends of the shaft (input/output is \theta_{ii} and meshing rotation is \theta_{ij}).
Parameters
| Newmark parameters | Value | Unit |
|---|---|---|
\gamma |
0.7 | |
\beta |
0.3 | |
| Gears | Value | |
Pinon 1 tooth number Z1 |
72 | |
Pinion 2 tooth number Z2 |
18 | |
Pinion modulus m_0 |
0.016 | m |
Pinion width b |
0.1 | m |
Internal diameter pinion 1 d_{i1} |
0.5 | m |
Internal diameter pinion 2 d_{i_2} |
0.3 | m |
Input speed N_1 |
17 | m/s |
Contact angle \alpha |
20 | degree |
Steel density \rho_{steel} |
7860 | kg/m3 |
| Masses | Value | Unit |
Input mass m_{in} |
54000 | kg |
Output mass m_{out} |
10 | kg |
| Stiffness | ||
Flexion x : k_x=k_{x_1}=k_{x_2} |
1e8 | N/m |
Flexion y : k_y=k_{y_1}=k_{y_2} |
1e8 | N/m |
Torsion z : k_\theta=k_{\theta_1}=k_{\theta_2} |
1e5 | N/m |
Mean mesh stiffness k_{mesh_{mean}} |
1e8 | N/rad |
| Damping | Value | Unit |
Mass related \eta |
5e-2 | N/(m/s)/(kg) |
Stiffness related \beta |
1e-6 | N/(m/s)/(N/rad) |
| Wind excitation | Value | Unit |
Blade radius r_{blade} |
6 | m |
Air density \rho_{air} |
1.225 | kg/m3 |
Wind speed (mean) v_{wind_{mean}} |
10 | m/s |
Fluctuating external torque T_{fluct} |
50 | N |
Fluctuating external frequency f_{fluct} |
6 | Hz |
Performance efficiency cp |
16/27 |
Energy equations
\frac{\partial Ep}{\partial q_i}+
\frac{d}{dt}\left(\frac{Ec}{\dot{q}i}\right)-
\frac{\partial Ec}{\partial q_i}+
\frac{\partial D}{\partial \dot{q}i}
=\frac{W}{\partial q_i}
\
Ec = \frac{1}{2}m_1\dot{x}1^2+
\frac{1}{2}m_1\dot{y}1^2+
\frac{1}{2}I{11}\dot{\theta}{11}^2+
\frac{1}{2}I{12}\dot{\theta}{12}^2+
\frac{1}{2}m_2\dot{x_2}^2+
\frac{1}{2}m_2\dot{y_2}^2+
\frac{1}{2}I_{22}\dot{\theta}{22}^2+
\frac{1}{2}I{21}\dot{\theta}{21}^2
\
Ep = \frac{1}{2}k{x_1}x_{1}^2+
\frac{1}{2}k_{y_1}y_{1}^2+
\frac{1}{2}k_{\theta_1}(\theta_{11}-\theta_{12})^2+
\frac{1}{2}k_{x_2}x_{2}^2+
\frac{1}{2}k_{y_2}y_{2}^2+
\frac{1}{2}k_{\theta_2}(\theta_{21}-\theta_{22})^2+
\frac{1}{2}K_e(t)\delta_1^2(t)
Equations of movement
Ep and Ec are separated between the two shaft (1 and 2)
\begin{cases}
&\frac{\partial Ep_1}{\partial x_1}+\frac{d}{dt}\left(\frac{\partial Ec_1}{\partial \dot{x}1}\right)=m_1\ddot{x}1+k{x1}x_1+sin(\alpha)K_e(t)\delta 1(t)=0 \
&\frac{\partial Ep_1}{\partial y_1}+\frac{d}{dt}\left(\frac{\partial Ec_1}{\partial \dot{y}1}\right)=m_1\ddot{y}1+k{y1}y_1+cos(\alpha)K_e(t)\delta 1(t)=0 \
&\frac{\partial Ep_1}{\partial \theta{11}}+\frac{d}{dt}\left(\frac{\partial Ec_1}{\partial \dot{\theta}{11}}\right)=I{11}\ddot{\theta}{11}+k_{\theta 1}(\theta_{11}-\theta_{12})=Cm \
&\frac{\partial Ep_1}{\partial \theta_{12}}+\frac{d}{dt}\left(\frac{\partial Ec_1}{\partial \dot{\theta}{12}}\right)=I{12}\ddot{\theta}{12}+k{\theta 1}(\theta_{12}-\theta_{11})+K_e(t)r_{b12}\delta_1(t)=0 \
&\frac{\partial Ep_1}{\partial x_2}+\frac{d}{dt}\left(\frac{\partial Ec_2}{\partial \dot{x}2}\right)=m_2\ddot{x}2+k{x2}x_2-sin(\alpha)K_e(t)\delta 1(t)=0 \
&\frac{\partial Ep_1}{\partial y_2}+\frac{d}{dt}\left(\frac{\partial Ec_2}{\partial \dot{y}2}\right)=m_2\ddot{y}2+k{y2}y_2-cos(\alpha)K_e(t)\delta 1(t)=0 \
&\frac{\partial Ep_1}{\partial \theta{22}}+\frac{d}{dt}\left(\frac{\partial Ec_2}{\partial \dot{\theta}{22}}\right)=I{22}\ddot{\theta}{22}-k_{\theta 2}(\theta_{21}-\theta_{22})=-Cr \
&\frac{\partial Ep_1}{\partial \theta_{21}}+\frac{d}{dt}\left(\frac{\partial Ec_2}{\partial \dot{\theta}{21}}\right)=I{21}\ddot{\theta}{21}-k{\theta 2}(\theta_{21}-\theta_{22})+K_e(t)r_{b21}\delta_1(t)=0
\end{cases}
With \delta_i(t)=(x_1-x_2)sin(\alpha)+(y_1-y_2)cos(\alpha)+\theta_{12}r_{b1}+\theta_{21}r_{b2} and r_{bij} the base radius of pinions, K_e is the varying meshing stiffness, mass m_i = \pi r^2_{i}\rho for both pinions, and I_{ij} = \frac{1}{2}m_{ij}r_{ij}^2 for pinions and input/output inertia.
Matrices
Equation to solve :
[M]{\ddot{q}}+C{\dot{q}}+([Ks]+[K(t)]){q}={F_0} \
{q}={x_1,y_1,\theta_{11},\theta_{12},x_2,y_2,\theta_{22},\theta_{21}}
Mass matrix :
M =
\begin{bmatrix}
m_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & m_1 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & I_{11} & 0 & 0 & 0 & 0 & 0 \
0 & 0 & 0 & I_{12} & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & m_2 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & m_2 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & I_{22} & 0 \
0 & 0 & 0 & 0 & 0 & 0 & 0 & I_{21} \
\end{bmatrix} \
With m_1 mass of shaft 1 and m_2 mass of shaft 2
Variable mesh stiffness matrix
K(t)= K_e(t)
\begin{bmatrix}
sin(\alpha)^2 & sin(\alpha).cos(\alpha) & 0 & r_1.sin(\alpha) & -sin(\alpha)^2 & -sin(\alpha).cos(\alpha) & 0 & r_2.sin(\alpha) \
sin(\alpha).cos(\alpha) & cos(\alpha)^2 & 0 & r_1.cos(\alpha) & -sin(\alpha).cos(\alpha) & -cos(\alpha)^2 & 0 & r_2.cos(\alpha) \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
r_1.sin(\alpha) & r_1.cos(\alpha) & 0 & r_1^2 & -r_1.sin(\alpha) & -r_1.cos(\alpha) & 0 & r_1.r_2 \
-sin(\alpha)^2 & -sin(\alpha).cos(\alpha) & 0 & -r_1.sin(\alpha) & sin(\alpha)^2 & sin(\alpha).cos(\alpha) & 0 & -r_2.sin(\alpha) \
-sin(\alpha).cos(\alpha) & -cos(\alpha)^2 & 0 & -r_1.cos(\alpha) & sin(\alpha).cos(\alpha) & cos(\alpha)^2 & 0 & -r_2.cos(\alpha) \
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
r_2.sin(\alpha) & r_2.cos(\alpha) & 0 & r_1.r_2 & -r_2.sin(\alpha) & -r_2.cos(\alpha) & 0 & r_2^2 \
\end{bmatrix}
Note
In matlab, cos and sin terms will be noted as
s_i, for visual simplification
Constant stiffness matrix
Ks=
\begin{bmatrix}
k_{x1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \
0 & k_{y1} & 0 & 0 & 0 & 0 & 0 & 0 \
0 & 0 & k_{\theta1}& -k_{\theta1} & 0 & 0 & 0 & 0 \
0 & 0 & -k_{\theta1} & k_{\theta1} & 0 & 0 & 0 & 0 \
0 & 0 & 0 & 0 & k_{x1} & 0 & 0 & 0 \
0 & 0 & 0 & 0 & 0 & k_{y1} & 0 & 0 \
0 & 0 & 0 & 0 & 0 & 0 & k_{\theta2} & -k_{\theta2} \
0 & 0 & 0 & 0 & 0 & 0 & -k_{\theta2} & k_{\theta2} \
\end{bmatrix}
Damping matrix
C = \eta M+\beta K_{mean}
Forces
F =
\begin{bmatrix}
0 \ 0 \ T_{in} \ 0 \ 0\ 0\ T_{out}\ 0
\end{bmatrix}
Mesh stiffness model
The mesh stiffness is modelized by a trapezoïdal signal between a low value (=k_{mesh_{mean}} for one tooth in contact) and its double value as a high value. Length of both phase is given by T_{high}=(\varepsilon-1)T_{mesh} and T_{low}=(2-\varepsilon)T_{mesh}. 10% of the start and of the end of the “high” phase is a linear slope from low value to high value.
Note
T_{mesh}=1/f_{mesh}=1/(Z1*N1/60)
\varepsilonis an approximation using the following formula :c_1=\frac{1+\frac{1}{Z1}}{\left(\pi cos(\alpha)\left(\frac{sin(\alpha)}{2}+\sqrt{\frac{sin(\alpha)^2}{4}+\frac{1}{Z1^2}+\frac{1}{Z1}}\right)\right)} \ c_1=\frac{1+\frac{1}{Z2}}{\left(\pi cos(\alpha)\left(\frac{sin(\alpha)}{2}+\sqrt{\frac{sin(\alpha)^2}{4}+\frac{1}{Z2^2}+\frac{1}{Z2}}\right)\right)} \ c_{12}=c_1+c_2\approx \epsilon
Model
Time : time is sampled on the T_{mesh}, with fs the number of points for each period, and n_{periods} the number of periods calculated. The number of points is fs\times n_{periods}.
Speed : rotational speed is initialized at its “stable” value.
if no automatic implementation, it has to be done “by hand” by executing the program once, get the speeds after the transient zone, and put them by hand in “initialization” source code (
func.m)
Initialization : everything is set at 0 for initialization (except speed)
Model : a Newmark scheme is used to compute the
TO FIX
Caution
...
Machine specifications
Example of computational time : 44sec for fs=2500 / n_periods = 10
-`
.o+`
`ooo/ OS: Arch Linux x86_64
`+oooo: Host: 20EV000UFR ThinkPad E560
`+oooooo: Kernel: 6.19.8-arch1-1
-+oooooo+: Shell: zsh 5.9
`/:-:++oooo+: Terminal: kitty
`/++++/+++++++: CPU: Intel i5-6200U (4) @ 2.300GHz
`/++++++++++++++: GPU: Intel Skylake-U GT2 [HD Graphics 520]
`/+++ooooooooooooo/` Memory: 1703MiB / 3792MiB
./ooosssso++osssssso+`
.oossssso-````/ossssss+`
-osssssso. :ssssssso.
:osssssss/ osssso+++.
/ossssssss/ +ssssooo/-
`/ossssso+/:- -:/+osssso+-
`+sso+:-` `.-/+oso:
`++:. `-/+/
.` `/