18 KiB
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
cd wind_turbine_gearbox/octave/
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 objects)
- 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)
- Display
- time response on a2
- speed response
- x movement (vs y movement ?)
- RMS STD Kurtosis
- Printing pdf and svg
- Convergence verification
- Display effect of defects
- Verify convergence for other input
- Preparation Monte Carlo Parallel calculation
- Decide a number of samples per mesh period
- Random models implementation
- 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}} |
1e8 | 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)\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 \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 (
+\inftyis approximated as1000) :
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
\alphadepending on a variation ofa, ingearbox_geometry. It does only change the values of contact ratio\varepsilonand 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:
`++:. `-/+/
.` `/