Pu_add estimation of alpha for more accurate radius of race at point of contact

@BearImp/calcCap.m

@@ -214,36 +214,36 @@ end
 
 
 function leander_parallel(indices)
-    C.R_hi = @(v) G.R_bi.*sqrt(1-(tempR_RE/G.R_bi.*cos(v)).^2);
-    C.R_ha = @(v) G.R_ba.*sqrt(1-(tempR_RE/G.R_ba.*cos(v)).^2);
+    C.R_hi = @(v) G.R_ri.*sqrt(1-(tempR_RE/G.R_ri.*cos(v)).^2);
+    C.R_ha = @(v) G.R_ra.*sqrt(1-(tempR_RE/G.R_ra.*cos(v)).^2);
     
     for ii = indices
-        C.I_i = @(v,p) abs(tempR_RE^2.*sin(v))./(-(C.R_hi(v)-tempR_li).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_hi(v)+tempR_li)).^2)+tempR_RE.*(1+sin(v).*cos(p))-tempR_li+G.R_bi+B.s(ii));
-        C.I_a = @(v,p) abs(tempR_RE^2.*sin(v))./( (C.R_ha(v)+tempR_la).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_ha(v)+tempR_la)).^2)+tempR_RE.*(1-sin(v).*cos(p))-tempR_la-G.R_ba+B.s(ii));
+        C.I_i = @(v,p) abs(tempR_RE^2.*sin(v))./(-(C.R_hi(v)-tempR_li).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_hi(v)+tempR_li)).^2)+tempR_RE.*(1+sin(v).*cos(p))-tempR_li+G.R_ri+B.s(ii));
+        C.I_a = @(v,p) abs(tempR_RE^2.*sin(v))./( (C.R_ha(v)+tempR_la).*sqrt(1-(tempR_RE.*sin(p).*sin(v)./(C.R_ha(v)+tempR_la)).^2)+tempR_RE.*(1-sin(v).*cos(p))-tempR_la-G.R_ra+B.s(ii));
         C.C_out(1,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(C.I_i,C.phi_1i(ii)+pi/2,C.phi_2i+pi/2, pi,1.5*pi);
         C.C_out(2,ii) = 4 * obj.epsilon_0 * S.epsilon_Oel * integral2(C.I_a,C.phi_1a(ii)+pi/2,C.phi_2a+pi/2,  0,0.5*pi);
     end
 end
 
 function leander_radial(indices)
-    C.c_i(1,indices) =  tempR_li - G.R_bi - tempR_RE - temp_s(1,indices);          % TODO: B.s an der Stelle sinnvoll? ###
-    C.c_a(1,indices) = -tempR_la - G.R_ba + tempR_RE + temp_s(2,indices);          % TODO: B.s an der Stelle sinnvoll? ###
+    C.c_i(1,indices) =  tempR_li - G.R_ri - tempR_RE - temp_s(1,indices);          % TODO: B.s an der Stelle sinnvoll? ###
+    C.c_a(1,indices) = -tempR_la - G.R_ra + tempR_RE + temp_s(2,indices);          % TODO: B.s an der Stelle sinnvoll? ###
     
     for ii = indices
-        C.I_i = @(g,t) abs(tempR_li*(tempR_li.*sin(g)+G.R_bi))./(sqrt(tempR_li.^2+G.R_bi.^2+C.c_i(ii').^2+2*G.R_bi*tempR_li.*sin(g)+2*C.c_i(ii')*(tempR_li.*sin(g)+G.R_bi).*cos(t))-tempR_RE);
-        C.I_a = @(g,t) abs(tempR_la*(tempR_la.*sin(g)+G.R_ba))./(sqrt(tempR_la.^2+G.R_ba.^2+C.c_a(ii').^2+2*G.R_ba*tempR_la.*sin(g)+2*C.c_a(ii')*(tempR_la.*sin(g)+G.R_ba).*cos(t))-tempR_RE);
+        C.I_i = @(g,t) abs(tempR_li*(tempR_li.*sin(g)+G.R_ri))./(sqrt(tempR_li.^2+G.R_ri.^2+C.c_i(ii').^2+2*G.R_ri*tempR_li.*sin(g)+2*C.c_i(ii')*(tempR_li.*sin(g)+G.R_ri).*cos(t))-tempR_RE);
+        C.I_a = @(g,t) abs(tempR_la*(tempR_la.*sin(g)+G.R_ra))./(sqrt(tempR_la.^2+G.R_ra.^2+C.c_a(ii').^2+2*G.R_ra*tempR_la.*sin(g)+2*C.c_a(ii')*(tempR_la.*sin(g)+G.R_ra).*cos(t))-tempR_RE);
         C.C_out(1,ii,posBall_conductive) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(C.I_i,-pi/2 + C.g_1i(ii,posBall_conductive),   -C.g_2i, 0, pi/9);
         C.C_out(2,ii,posBall_conductive) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(C.I_a, pi/2 + C.g_1a(ii,posBall_conductive), pi-C.g_2a, 0, pi/9);
     end 
 end
 
 function leander_steffen(indices)
-    C.dA_i = @(p) tempR_li * (G.R_bi + tempR_li*cos(p));
-    C.dA_a = @(p) tempR_la * (G.R_ba + tempR_la*cos(p));
-    C.x_kmi = G.R_bi + C.DeltaR_i;
-    C.x_kma = G.R_ba - C.DeltaR_a;
-    C.h_i = @(p,t,ii) sqrt(G.R_bi^2 + tempR_li^2 + C.x_kmi(ii)^2 + 2*G.R_bi*tempR_li*cos(p) - 2*C.x_kmi(ii)*cos(t).*(G.R_bi+tempR_li*cos(p)))  -  tempR_RE;
-    C.h_a = @(p,t,ii) sqrt(G.R_ba^2 + tempR_la^2 + C.x_kma(ii)^2 + 2*G.R_ba*tempR_la*cos(p) - 2*C.x_kma(ii)*cos(t).*(G.R_ba+tempR_la*cos(p)))  -  tempR_RE;
+    C.dA_i = @(p) tempR_li * (G.R_ri + tempR_li*cos(p));
+    C.dA_a = @(p) tempR_la * (G.R_ra + tempR_la*cos(p));
+    C.x_kmi = G.R_ri + C.DeltaR_i;
+    C.x_kma = G.R_ra - C.DeltaR_a;
+    C.h_i = @(p,t,ii) sqrt(G.R_ri^2 + tempR_li^2 + C.x_kmi(ii)^2 + 2*G.R_ri*tempR_li*cos(p) - 2*C.x_kmi(ii)*cos(t).*(G.R_ri+tempR_li*cos(p)))  -  tempR_RE;
+    C.h_a = @(p,t,ii) sqrt(G.R_ra^2 + tempR_la^2 + C.x_kma(ii)^2 + 2*G.R_ra*tempR_la*cos(p) - 2*C.x_kma(ii)*cos(t).*(G.R_ra+tempR_la*cos(p)))  -  tempR_RE;
     for ii = indices
         x1=double(C.g_1i(ii));
         x2=double(C.g_2i);
@@ -283,8 +283,8 @@ function runSemianalytisch3D(indices)
     
     if AddOn.Parallel_Computing_Toolbox
 
-        tmp_semi_i = @(runSemi) semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_bi, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
-        tmp_semi_a = @(runSemi) semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ba, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
+        tmp_semi_i = @(runSemi) semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_ri, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
+        tmp_semi_a = @(runSemi) semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ra, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
 
         parfor runSemi = indices
             C_semi_i(runSemi) = feval(tmp_semi_i, runSemi); %da parfor loop nicht auf nested functions (semianalytisch3D) direkt zugreifen kann
@@ -294,8 +294,8 @@ function runSemianalytisch3D(indices)
     else
 
         for runSemi = 1 : numel(temp_s)/2
-           C_semi_i(runSemi) = semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_bi, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
-           C_semi_a(runSemi) = semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ba, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
+           C_semi_i(runSemi) = semianalytisch3D(temp_s(1,runSemi), tempR_RE, tempR_li, G.B_i, L.B, -G.R_ri, S.epsilon_Oel,temp_a(1,runSemi),temp_b(1,runSemi),110,110);
+           C_semi_a(runSemi) = semianalytisch3D(temp_s(2,runSemi), tempR_RE, tempR_la, G.B_a, L.B,  G.R_ra, S.epsilon_Oel,temp_a(2,runSemi),temp_b(2,runSemi),110,110);
         end
     end
 
@@ -393,7 +393,6 @@ function C_Zusatz=semianalytisch3D(s,R_WK,R_R,B_R,B,R_L,epsilon_r,varargin)
     elseif  any(r2small)
 %         fprintf('%i Entries of r have been removed \n',sum(r2small))
         r(r2small)=inf;
-    
     end
     
     integrand_t = R_WK^2./(r-R_WK).*cos(T_t);

@BearImp/calcGeo.m

-function calcGeo(obj)
+function calcGeo(obj,options)
 %2.3 LagerGeometrie G = G(L,R,T_Oil)
 %Berechnet verschiedene Größen der Lagergeometrie wie z.B.
 %Berührungswinkel, Rollradien, Elliptizitätsparameter der Hertz'schen
@@ -7,31 +7,36 @@ function calcGeo(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler, Julius van der Kuip
 
+arguments
+    obj
+    options.Delta_Tia = 5 % Temperaturerhöhung des Innenrings gegenüber des Außenrings
+    options.alpha = []
+end
 %% Parameter prüfen
 assert(obj.up2date.L,'Lager nicht gesetzt')
 assert(obj.up2date.R,'Lagerspiel noch nicht berechnet')
-L = obj.L; R = obj.R; T_Oil = obj.T_Oil;
-G = struct; 
+L = obj.L; R = obj.R; T_Oil = obj.T_Oil; F_a = obj.F_a;
+G = struct; G.method = obj.method.G;
 
 %% Berechnung
-G.alpha = 0; 
-
+G.Delta_Tia = options.Delta_Tia;                 % Temperaturerhöhung des Innenrings gegenüber des Außenrings
 G.Delta_Ta_norm = T_Oil - 20;                    % Temperaturerhöhung der Welle bezogen auf die Normtemperatur von 20°C (Innenring heißer als Außenring)
-G.Delta_Ti_norm = G.Delta_Ta_norm + 5; %Temperaturerhöhung des Gehäuses bezogen auf die Normtemperatur von 20°C #######Idee: 5°C durch festlegbare Variable ersetzen
+G.Delta_Ti_norm = G.Delta_Ta_norm + G.Delta_Tia; % Temperaturerhöhung des Gehäuses bezogen auf die Normtemperatur von 20°C
+G.Delta_Tm_norm = (G.Delta_Ta_norm + G.Delta_Ti_norm)/2;
 
 G.E_red(1) = 2  ./  ((1-L.RE_A.nu.^2)./L.RE_A.E + (1-L.ring.nu.^2)/L.ring.E);
-G.D_RE(1) = L.D_RE * (1 + L.RE_A.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
+G.D_RE(1) = L.D_RE * (1 + L.RE_A.alpha_t * G.Delta_Tm_norm);
 
-if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnet
+if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnen
     G.E_red(2) = 2  ./  ((1-L.RE_B.nu.^2)./L.RE_B.E + (1-L.ring.nu.^2)/L.ring.E);
-    G.D_RE(2) = L.D_RE * (1 + L.RE_B.alpha_t * ((G.Delta_Ti_norm + G.Delta_Ta_norm) / 2));
+    G.D_RE(2) = L.D_RE * (1 + L.RE_B.alpha_t * G.Delta_Tm_norm);
 end
 
 G.D = (L.D - R.e_ar) * (1 + L.ring.alpha_t * G.Delta_Ta_norm);
 G.d = (L.d + R.e_ir) * (1 + L.ring.alpha_t * G.Delta_Ti_norm);
+G.R_RE = G.D_RE/2;
 
 G.d_m = (G.d + G.D)/2;
-G.R_RE = G.D_RE/2;
 
 b_li = (L.D-L.d-2*L.D_RE-L.c)/2; %Breite von Innenring (bezogen auf Durchmesser)
 b_la = b_li; 
@@ -46,27 +51,57 @@ G.B_a=L.B_a*(1+L.ring.alpha_t*G.Delta_Ta_norm);
 G.R_li = L.f_i*G.D_RE; %Laufbahnradius des Innenrings (Dowson S.48) %der im Lagerquerschnitt gemessene Laufbahnrillenradius des Außenrings
 G.R_la = L.f_a*G.D_RE; %Laufbahnradius des Außenrings
 
+G.b = L.f_i + L.f_a - 1; % Harris B
+if ~isempty(options.alpha)
+    G.alpha = options.alpha;
+elseif ~F_a || G.method.alphaForRaceRadius == 'zero'
+    G.alpha = 0;
+elseif G.method.alphaForRaceRadius == 'estimate'
+    if G.P_d > 0
+        G.alpha_0 = acos(1 - G.P_d/(2*G.D_RE*G.b));
+    else
+        G.alpha_0 = 0;
+    end
+    G.k_ax_est = interp1([0 0.02 0.08 0.14 0.2],[0 50   250  500  760]*1e3,G.b,'spline')*6894.757; % Harris.2001 Fig. 7.5, p. 247
+    alpha = pi/3; % Stable start value
+    numberOfIterations = 1;
+    while true
+        cosAlphaRatioTerm = (cos(G.alpha_0)/cos(alpha) - 1);
+        alpha_next = alpha + (            F_a/(L.Z*L.D_RE^2*G.k_ax_est) - sin(alpha)*cosAlphaRatioTerm^1.5               )  /  ... Harris.2001 (7.34), p. 247
+                             ( cos(alpha)*cosAlphaRatioTerm^1.5  +  1.5*tan(alpha)^2*cosAlphaRatioTerm^0.5*cos(G.alpha_0) );
+        if abs(alpha-alpha_next) <= 10*eps(alpha)
+            break
+        elseif numberOfIterations >= 1000
+            warning('Estimation of mounted contact angle did not converge')
+            break
+        else
+            alpha = alpha_next;
+            numberOfIterations = numberOfIterations + 1;
+        end
+    end
+    G.alpha = alpha_next;
+end
+
 G.R_bi = (G.d + G.b_li) / 2; %Minimaler innerer Rillenradius
-G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius
+G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius                 %  unabh. von alpha
 
-G.sigmarho_i = (2./G.R_RE + 1./G.R_bi - 1./G.R_li);
-G.sigmarho_a = (2./G.R_RE - 1./G.R_ba - 1./G.R_la);
+G.R_ri = G.R_bi + G.R_RE*(1-cos(G.alpha)); % Race radius at point of contact
+G.R_ra = G.R_ba - G.R_RE*(1-cos(G.alpha));
 
-G.costau_i = (1./G.R_bi + 1./G.R_li) ./ G.sigmarho_i;
-G.costau_a = (-1./G.R_ba + 1./G.R_la) ./ G.sigmarho_a;
+G.sigmarho_i = (2./G.R_RE + 1./G.R_ri - 1./G.R_li);
+G.sigmarho_a = (2./G.R_RE - 1./G.R_ra - 1./G.R_la);                         % abh. von alpha
+G.R_si = 1 ./ G.sigmarho_i;                                                                     % in calcFilm
+G.R_sa = 1 ./ G.sigmarho_a;                                                 % abh. von alpha    % in calcFilm
 
-G.R_xi = 1 ./ ( 1./G.R_RE + 1./G.R_bi );
-G.R_xa = 1 ./ ( 1./G.R_RE - 1./G.R_ba );
-G.R_x  = [G.R_xi;G.R_xa];
+G.R_xi = 1 ./ ( 1./G.R_RE + 1./G.R_ri );
+G.R_xa = 1 ./ ( 1./G.R_RE - 1./G.R_ra );
+G.R_x  = [G.R_xi;G.R_xa];                                                   % abh. von alpha    % in calcFilm & calcCap
 
 G.R_yi = 1 ./ ( 1./G.R_RE - 1./G.R_li );
 G.R_ya = 1 ./ ( 1./G.R_RE - 1./G.R_la );
 G.R_y  = [G.R_yi;G.R_ya];
 
-G.R_si = 1 ./ G.sigmarho_i;
-G.R_sa = 1 ./ G.sigmarho_a;
-
-G.gamma = G.D_RE  *cos(G.alpha) ./ G.d_m;
+G.gamma = G.D_RE  *cos(G.alpha) ./ G.d_m;                                   % abh. von alpha   % in calcFilm
 
 integralF = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^-0.5;
 integralE = @(phi,k) (1 - (1-1/k.^2).*sin(phi).^2).^0.5;
@@ -75,23 +110,23 @@ fraktF = @(k) integral( @(phi) integralF(phi,k),0,pi/2 );
 fraktE = @(k) integral( @(phi) integralE(phi,k),0,pi/2 );
 
 G.F_i = (( G.gamma/(1-G.gamma)+G.D_RE/(2*G.R_li))/(2+G.gamma/(1-G.gamma)-G.D_RE/(2*G.R_li)));
-G.F_a = ((-G.gamma/(G.gamma+1)+G.D_RE/(2*G.R_la))/(2-G.gamma/(1+G.gamma)-G.D_RE/(2*G.R_la)));
+G.F_a = ((-G.gamma/(1+G.gamma)+G.D_RE/(2*G.R_la))/(2-G.gamma/(1+G.gamma)-G.D_RE/(2*G.R_la)));     % abh. von alpha
 
 k_i_fkt = @(k_i) (1-2/(k_i^(2)-1)*((fraktF(k_i)/fraktE(k_i))-1)-G.F_i);
-k_a_fkt = @(k_a) (1-2/(k_a^(2)-1)*((fraktF(k_a)/fraktE(k_a))-1)-G.F_a);
+k_a_fkt = @(k_a) (1-2/(k_a^(2)-1)*((fraktF(k_a)/fraktE(k_a))-1)-G.F_a);     % abh. von alpha
 
-G.k_i = fzero(k_i_fkt,5);
-G.k_a = fzero(k_a_fkt,5);
-G.k = [G.k_i;G.k_a];
+G.k_i = fzero(k_i_fkt,5);                                                                       % in calcFilm
+G.k_a = fzero(k_a_fkt,5);                                                                       % in calcFilm
+G.k = [G.k_i;G.k_a];                                                        % abh. von alpha    % in calcFilm & calcCap
 
 G.fraktF_i = fraktF(G.k_i);
 G.fraktF_a = fraktF(G.k_a);
-G.fraktE_i = fraktE(G.k_i);
-G.fraktE_a = fraktE(G.k_a);
+G.fraktE_i = fraktE(G.k_i);                                                                     % in calcFilm
+G.fraktE_a = fraktE(G.k_a);                                                 % abh. von alpha    % in calcFilm
 
-G.K_pi = pi.*G.k_i.*G.E_red.*(G.fraktE_i./(4.5*G.sigmarho_i.*G.fraktF_i.^3)).^0.5;
-G.K_pa = pi.*G.k_a.*G.E_red.*(G.fraktE_a./(4.5*G.sigmarho_a.*G.fraktF_a.^3)).^0.5;
-G.K_p= 1./(1./G.K_pi.^(2/3)+1./G.K_pa.^(2/3)).^1.5;
+G.K_pi = pi.*G.k_i.*G.E_red.*(G.fraktE_i./(4.5*G.sigmarho_i.*G.fraktF_i.^3)).^0.5;              % in calcLoad
+G.K_pa = pi.*G.k_a.*G.E_red.*(G.fraktE_a./(4.5*G.sigmarho_a.*G.fraktF_a.^3)).^0.5;              % in calcLoad
+G.K_p= 1./(1./G.K_pi.^(2/3)+1./G.K_pa.^(2/3)).^1.5;                         % abh. von alpha    % in calcLoad
 
 G.Deltapsi = 2*pi./L.Z;
 G.psi_ = (0:L.Z-1) .* G.Deltapsi;

@BearImp/calcLoad.m

@@ -206,7 +206,7 @@ end
 
 function Q_total_y_perball_cont=totalReactionForcePerBallOnlyY(G,IndBall,delta_s,psi_run)
    delta_nury = (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))/cos(atan(sin(psi_run(IndBall))*delta_s(1)/ ...
-                (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))))-G.R_RE-G.R_bi - G.D_RE; %- G.D_RE richtig?
+                (G.R_ba-G.R_RE-cos(psi_run(IndBall))*delta_s(1))))-G.R_RE-G.R_bi - G.D_RE;
    Q_total_y_perball_cont = abs((delta_nury<0)*(delta_nury))^(3/2)*G.K_p;
 end
 

BearImpOptions.m

@@ -28,6 +28,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
 
     methods (Static,Hidden) 
         function possibleOption = Possible
+            possibleOption.G.alphaForRaceRadius = {'zero','estimate'};
             possibleOption.T = {'linear','Vogel'};
             possibleOption.B = {'static','dynamic'};
             possibleOption.H.filmThicknessFormula = {'Hamrock/Dowson','Moes'};
@@ -43,6 +44,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
 
         function defaultOption = Default
         % Gibt ein Rechenmethoden-Struct mit den Default-Methoden aus
+            defaultOption.G.alphaForRaceRadius = 'estimate';
             defaultOption.T = 'Vogel';
             defaultOption.B = 'dynamic';
             defaultOption.H.filmThicknessFormula = 'Hamrock/Dowson';