--- Author: Antoine Foucault-Castelli Date: 23/03/2026 Update: 31/03/2026 Bibliography: @karmiFiabiliteOptimisationSystemes @ghorbelEffectBrakeLocation2020a @zhengAnalyticalApproachMesh2022 Version: 5.0 --- > [!BUG] > Some equations has been put in latex block code as they didnt show properly on gitlab # Gearbox simulation ## Structure ``` %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -- main main.m >> used to launch all run_*.m at once main_parameters.m >> charge model's parameters -- runners run_main.m >> PURPOSE: launch every class_func functions in the right order to model the system run_display.m >> PURPOSE: launch every class_displayer functions run_print.m >> PURPOSE: print the figures made by run_display >> OUTPUT: .svg -- classes class_display.m >> displays the inputs and outputs of the model [col1,col2,col3,col4,col5,col6,col7,col8,lw1,lw2] = charge_aff() [figname,namex,namey,namet,ylabelx,ylabelt,ylabelfft,ylabelfftT] = type_choice(type) listOfFig = disp_x(t,x,type,listOfFig) listOfFig = excitations_ext(t,T_aero,k_var,k_mean,listOfFig) listOfFig = spectral_analysis(x,tf,dt,maxf,type,listOfFig) % TODO plot_fep(t,Fep,i) listOfFig = plotThetaDiff(t,x,listOfFig) class_func.m >> model functions (explicit cf. comments in code) [t, dt, tf, t_mesh, f_mesh, fs] = time_sampling(fs_div,n_periods,Z1,N1) [Jin, J1, J2, Jout, J_all] = inertia(m_all,r_blade,di1,di2) [c12, db1, db2] = gearbox_geometry(alpha0, Z1, Z2, m0, a) [k_var, k_mean] = mesh_stiffness(epsilon, t, t_mesh, k_mesh_mean) [T_aero,T_out] = wind_excitation(t,rho_air,r_blade,v_wind,T_fluct_in,f_ext,omega_in,cp,ratio) [M] = mass_matrix(m_all, J_all) [K_var, K_const, K_mean] = stiffness_matrix(db1, db2, alpha0, kx, ky, ktheta, k_mean) [K] = stiffness_update(K_const,K_var,k) [C] = damping_matrix(M,K_mean,damp_eta,damp_beta) [F] = forces_matrix(T_in,T_out) [x,v,a] = initialization(N1,t,ratio,cp) [x,v,a] = newmark(M,C,K_const,K_var,K_mean,k_var,dt,t,T_in,T_out,x,v,a,omega_in,gamma,beta,F_ep) [F_ep] = error_profile_excitation(e12,f_mesh,t,alpha0,db1,db2,k_var) [damp_eta, damp_beta] = damping_coeff(M,K,xi1,xiN) -- standalone _convergence.m >> purpose: test the numerical convergence of the model >> output: 3 graphs + print/_convergence/*.svg _gear-mesh-var.m >> purpose: shows the variation of gear mesh with central distance changes >> output: 1 graph + print/_gear-mesh-var/*.svg _rms-analysis.m >> purpose: analysis of RMS, STD and kurtosis for many setups >> output: data/rms.svg ``` ## How does it work ? **For Matlab under windows :** classic matlab > [!WARNING] > NOT TESTED UNDER MATLAB > remove all occurences of `graphics_toolkit("gnuplot")` in the code which is only useful for Octave visualization **For Octave/Linux :** *May need a few package installation (follow octave recommandations after first(s) run* ```baswind_turbine_gearbox/octave/h cd wind_turbine_gearbox/octave/ octave >>> run main.m >>> ``` ## Informations - [x] Verified code with base code and paper - [x] Implemented wind excitation - [x] Implemented Ghorbel parameters - [~] Add varying output torque depending on number of tooth in contact ? - [~] add brake ? - [~] Variation of C ? - [~] ODE45 ? - [x] Displays from last version (in class objects) - [x] New mesh stiffness variation model (trapezoidal) - [x] Print in pdf - [x] Errors/defects implementation - [x] profile error - [x] assembly defect (change on gear stiffness) - [x] README - [x] Maths - [~] Images - [x] Source(s) - [x] Display - [x] time response on a2 - [x] speed response - [x] x movement (vs y movement ?) - [x] RMS STD Kurtosis - [x] Printing pdf and svg - [x] Convergence verification - [x] Display effect of defects - [ ] Verify convergence for other input - [ ] Preparation Monte Carlo Parallel calculation - [x] Decide a number of samples per mesh period - [ ] Random models implementation - [x] Implemented easy disp and print - [ ] Add a "peak finder" for FFT ## 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 | | Rotor radius $r_{rot}$ | 0.5 | m | | Input theoretical speed $N_1$ | 17 | m/s | | Contact angle $\alpha$ | 20 | deg | | Steel density $\rho_{steel}$ | 7860 | kg/m3 | | **Masses** | **Value** | **Unit** | | Input mass $m_{in}$ (hub + blades) | 54000 | kg | | Output mass $m_{out}$ (rotor) | 10000 | 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}$ | 1e6 | N/m | | Mean mesh stiffness $k_{mesh_{mean}}$ | 2.8e8 | N/rad | | **Damping** | **Value** | **Unit** | | Mass damping ratio | 0.05 | | | Stiffness damping ratio | 0.01 | | | **Wind excitation** | **Value** | **Unit** | | Blade radius $r_{blade}$ | 6 | m | | Air density $\rho_{air}$ | 1.225 | kg/m3 | | Wind speed (mean) $v_{wind_{mean}}$ | 37.5 | m/s | | Fluctuating external torque $T_{fluct}$ | 50 | N | | Fluctuating external frequency $f_{fluct}$ | 6 | Hz | | Performance efficiency $cp$ | 16/27 | | | **Errors** | **Value** | **Unit** | | Profile error (approx) | 1e-5 | m | | Center distance error (approx) | 1e-2 | m | ### Energy equations $$ \left[\frac{\partial E_p}{\partial q_i}\right] + \left[\frac{d}{dt}\left(\frac{\partial E_k}{\partial\dot{q}_i}\right)\right] + \frac{\partial C}{\partial\dot{q}_i} = \frac{\partial W}{\partial q_i} \label{eq:reformlagrange} \\ E_k = \frac{1}{2} \left( m_1{\dot{x}_1}^2+ m_1{\dot{y}_1}^2+ I_{11}{\dot{\theta}_{11}}^2+ I_{12}{\dot{\theta}_{12}}^2+ m_2{\dot{x}_2}^2+ m_2{\dot{y}_2}^2+ I_{22}{\dot{\theta}_{22}}^2+ I_{21}{\dot{\theta}_{21}}^2 \right) \\ E_p = \frac{1}{2} \left( k_{x_1}{x_{1}}^2+ k_{y_1}{y_{1}}^2+ k_{\theta_1}(\theta_{11}-\theta_{12})^2+ k_{x_2}{x_{2}}^2+ k_{y_2}{y_{2}}^2+ k_{\theta_2}(\theta_{21}-\theta_{22})^2+ K_e(t){\delta}^2 \right) $$ ### Equations of movement The equations of movement are given by putting $E_k$ and $E_p$ into the Lagrangian formula. $$ \begin{cases} \frac{\partial E_p}{\partial x_1} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{x}_1} \right) + \frac{\partial C}{\partial \dot{x}_1} = m_1\ddot{x}_1 + k_{x_1}x_1 + \sin(\alpha)K_e(t)\delta + C_{x_1}\dot{x}_1 = 0 \\ \frac{\partial E_p}{\partial y_1} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{y}_1} \right) + \frac{\partial C}{\partial \dot{y}_1} = m_1\ddot{y}_1 + k_{y_1}y_1 + \cos(\alpha)K_e(t)\delta + C_{y_1}\dot{y}_1 = 0 \\ \frac{\partial E_p}{\partial \theta_{11}} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{\theta}_{11}} \right) + \frac{\partial C}{\partial \dot{\theta}_{11}} = I_{11}\ddot{\theta}_{11} + k_{\theta_{11}}(\theta_{11}-\theta_{12}) + C_{\theta_{11}}\dot{\theta}_{11} = T_{aero} \\ \frac{\partial E_p}{\partial \theta_{12}} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{\theta}_{12}} \right) + \frac{\partial C}{\partial \dot{\theta}_{12}} = I_{12}\ddot{\theta}_{12} - k_{\theta_{11}}(\theta_{11}-\theta_{12}) + {r_b}_1K_e(t)\delta + C_{\theta_{12}}\dot{\theta}_{12} = 0 \\ \frac{\partial E_p}{\partial x_2} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{x_2}} \right) + \frac{\partial C}{\partial \dot{x_2}} = m_2\ddot{x}_2 + k_{x_2}x_2 - \sin(\alpha)K_e(t)\delta + C_{x_2}\dot{x}_2 = 0 \\ \frac{\partial E_p}{\partial y_2} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{y_2}} \right) + \frac{\partial C}{\partial \dot{y_2}} = m_2\ddot{y}_2 + k_{y_2}y_2 - \cos(\alpha)K_e(t)\delta + C_{y_2}\dot{y}_2 = 0 \\ \frac{\partial E_p}{\partial \theta_{22}} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{\theta}_{22}} \right) + \frac{\partial C}{\partial \dot{\theta}_{22}} = I_{22}\ddot{\theta}_{22} - k_{\theta_{22}}(\theta_{21}-\theta_{22}) + C_{x_1}\dot{\theta}_{22} = T_{gen} \\ \frac{\partial E_p}{\partial \theta_{21}} + \frac{d}{dt}\left( \frac{\partial E_k}{\partial \dot{\theta}_{21}} \right) + \frac{\partial C}{\partial \dot{\theta}_{21}} = I_{21}\ddot{\theta}_{21} + k_{\theta_{21}}(\theta_{21}-\theta_{22}) - {r_b}_2K_e(t)\delta + C_{\theta_{21}}\dot{\theta}_{21} = 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 :** $$ \left[M\right]\{\ddot{q}\} + \left[C\right] \{\dot{q}\}+ \left(\left[K_{var}\right]+\left[K_{cst}\right]\right) \{q\} = \left\{F\right\} $$ **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_{var} = K_e(t) \begin{bmatrix} s^2 & sc & 0 & sr_1 & -s^2 & -sc & 0 & -sr_2 \\ cs & c^2 & 0 & cr_1 & -cs & -c^2 & 0 & -cr_2 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ r_1s & r_1c & 0 & r_1^2 & -r_1s & -r_1c & 0 & -r_1r_2 \\ -s^2 & -sc & 0 & -sr_1 & s^2 & sc & 0 & sr_2 \\ -cs & -c^2 & 0 & -cr_1 & cs & c^2 & 0 & cr_2 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -r_2s & -r_2c & 0 & -r_2r_1 & r_2s & r_2c & 0 & r_2r_1 \end{bmatrix} $$ With - $s=sin(\alpha)$ - $c=cos(\alpha) $ - $r1=r_{b_1} $ - $r2=r_{b_2}$ **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** Damping matrix is : $C = \eta M+\beta K_{mean}$ **Forces** $$ F = \begin{bmatrix} 0 \\ 0 \\ T_{aero} \\ 0 \\ 0\\ 0\\ T_{gen}\\ 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)$ > $\varepsilon$ is 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 \varepsilon $$ ## Numerical model, wind model and errors **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`) **Wind :** cosine variation around mean torque value $$ T_{in} = \rho_{air}\pi(r_{blade}^2)(v_{wind}^3)*Cp)/(2\omega) \\ T_{aero}(t) = T_{in} + T_{fluct}\cos(2\pi f_{ext}t) $$ **Initialization :** everything is set at 0 for initialization (except speed) **Model :** a basic Newmark scheme is used to compute the model over time **Error implementation :** - error profile is implemented as ($+\infty$ is approximated as $1000$) : $$ ep(t) = e_{12}+\sum ^{+\infty}_{n=1}e_{12}sin(2n\pi f_m t) \\ \left\{F_{ep}\right\} = \frac{\partial \delta (t)}{\partial q_i}k(t)ep(t) \\ \frac{\partial \delta(t)}{\delta q_i} = \begin{bmatrix} sin(\alpha) \\ cos(\alpha) \\ 0 \\ rb_1 \\ -sin(\alpha) \\ -cos(\alpha) \\ 0 \\ rb_2 \end{bmatrix} $$ - center distance error is implemented as a transformation of $\alpha$ depending on a variation of $a$, in `gearbox_geometry`. It does only change the values of contact ratio $\varepsilon$ and then the profile of the mesh stiffness (length of phases and max/min values). $$ \alpha'=acos\left(\frac{rb_1+rb_2}{E+a}\right) $$ ## TO FIX > [!CAUTION] > Speed not right ## Machine specifications For the computational speed, here’s my machine’s spec. : ``` -` .o+` `ooo/ OS: Arch Linux x86_64 `+oooo: Host: 20EV000UFR ThinkPad E560 `+oooooo: Kernel: 6.19.10-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: `++:. `-/+/ .` `/ ```