
% Données du cahier des charges

C=20; omega=721*pi/30; omegamax=1442*pi/30;

% Propriétés des matériaux

kr=0.5; rocu=1.72e-8; qt=2.5; Br=1.045; dt=7850; ft=50;
Bc=0.05; da=7400; Bt=1.5; mua=1.05; dcu=8950; kfoi=0.95;
muo=4*pi*1e-7; dcr=7850; Pm=15; kfui=0.8;

% Propriétés des matériaux pour la thermique

alphaa=-5e-4; alphacu=3.8e-3; Text=50; h=10;

% Paramètres d'entrée fixés

e=0.8e-3; Udc=120; Bcr=1.2; rrs=1.11; p=6; Ne=3/2*p;

% Paramètres d'entrée à trouver

Bd=1.8; Be=0.75; Bcs=0.8; Lm=45e-3; Ds=189e-3; delta=3e6;

% Equations dimensionnantes

alpha=pi/p;
beta=pi/p;
alphai=alpha/5;
kv=omegamax/omega;
E=Udc/2/kv;
I=C*omega/2/E;
n=4*E/Be/Ds/Lm/omega;
li=Be/Bd*alphai*Ds/2;
ld=Be/Bd*alpha*Ds/2;
eb=(Be/Bd*(alpha/2-asin(ld/Ds))*Ds/2)+Ds/2*(1-cos(asin(ld/Ds)));
hcs=1/2*Bd/Bcs*ld;
Ba=Be*alpha*Ds/2/beta/(Ds/2+e)/rrs/kfui;
hcr=1/2*Ba/Bcr*beta*(Ds/2+e);
Senc=3/2*n*I/delta/kr;

% calcul discriminant pour hauteur dent

a=pi;
b=-(2*pi*(Ds/2-eb)-Ne*(li+ld));
c=Senc;
discr=b^2-4*a*c;

% calcul hauteur dent

hd=(-b-sqrt(abs(discr)))/2/a;
Dint=Ds-2*(eb+hd+hcs);
f=p*omega/2/pi;
hc=(eb/cos(alpha/2)-Ds/2*(1/cos(alpha/2)-1));
hi=(Ds/2*(1-cos(alphai/2))+hc*cos(alphai/2));
Rtb=((ld-li)/4+(Ds/2-eb-hd/2)*pi/2/Ne);
Lds=(Lm/kfoi+pi*Rtb);
Ltot=(Lm/kfoi+2*((Ds/2-eb-hd/2)*pi/Ne-li/2));
Mcs=dt*pi*hcs*(-hcs+2*(Ds/2-eb-hd))*Lm;
Mds=dt*Ne*((ld+li)*hd+Ds/2*(alphai*(hi+hc)/2+alpha*(eb+hc)/2))*Lm;
Pf=qt*((f/ft)^1.5)*(Mcs*((Bcs/Bt)^2)+Mds*((Bd/Bt)^2));

% système d'équations implicites

Ta=100; Tcu=110;
while abs(Ta-Tcu)/Ta>1e-6
  Ta=Tcu;
  ha=-mua*Be*e/(Ba-Br*(1+alphaa*Ta));
  Dext=Ds+2*(e+ha+hcr);
  Rph=rocu*(1+alphacu*Tcu)*n/2*Lds*delta/I;
  Pj=2*Rph*I^2;
  Sext=pi/2*(Dext^2)+pi*Dext*Ltot;
  Tcu=Text+(Pj+Pf+Pm)/h/Sext;
end

% fin système d'équations implicites

Imax=-4*p/n*(Bc/muo*beta/alpha*(1+2*e/Ds)*rrs*kfui*e+(Bc-Br*(1+alphaa*Ta))*ha/muo/mua);
Ma=da*p*beta*ha*(ha+2*(Ds/2+e))*Lm*rrs;
Mcr=dcr*pi*hcr*(hcr+2*(Ds/2+e+ha))*Lm*rrs;
Mcu=dcu*3/2*n*I/delta*Lds;
Mtot=Ma+Mcr+Mcs+Mds+Mcu;
eta=(C*omega-Pm)/(C*omega+Pj+Pf);
FencIrc=muo/16*(n^2)/Ne*hd*Lm/((Ds/2-eb-hd/2)*pi/Ne-(li+ld)/2);
FentIrc=muo*3/16*(n^2)/Ne*Lm/(e+ha)*alpha*Ds/2;
FbecIrc=muo*3/16*(n^2)/Ne*hc*Lm/(Ds/2-hc/2)/(pi/Ne-(alpha+alphai)/2);
Lph=3/2*FentIrc+2*(FencIrc+FbecIrc);
t1=-Lph/Rph*log(abs(1-3*Rph*I/Udc/(2-1/kv)));


