From 92e5bb027bc6effbaaff1ce022ef036dd1222fa7 Mon Sep 17 00:00:00 2001
From: sp89hili <sp89hili@pmd196.csi.tu-darmstadt.de>
Date: Wed, 22 Jan 2025 09:17:39 +0100
Subject: [PATCH] Pu_Add options T='DIN 51563', C.roughnessCorrection='Napel'
 and C.rhoRatio='Tait', allow to use set method for BearImpOptions to level 3,
 allow to use calcLub, calcGeo and calcClear with arrays, add BearImp.empty
 for array creation

---
 @BearImp/BearImp.m   | 15 +++++++++++++++
 @BearImp/calcCap.m   | 25 +++++++++++++++++++------
 @BearImp/calcClear.m |  5 +++++
 @BearImp/calcGeo.m   |  6 ++++++
 @BearImp/calcLub.m   | 25 +++++++++++++++++++------
 @BearImp/set.m       | 18 ++++++++++++++----
 BearImpOptions.m     |  6 ++++--
 7 files changed, 82 insertions(+), 18 deletions(-)

diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index 4119c65..6f8b6ef 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -422,4 +422,19 @@ classdef BearImp < handle & matlab.mixin.Copyable
             obj = s;
         end
     end
+    methods (Static)
+        function obj = empty(n,m)
+            % Creates a nxm array of empty BearImp objects
+            arguments
+                n {mustBeInteger, mustBeScalarOrEmpty}
+                m {mustBeScalarOrEmpty} = n
+            end
+            obj(n,m) = BearImp;
+            for ii = 1:n
+                for jj = 1:m
+                    obj(ii,jj) = BearImp;
+                end
+            end
+        end
+    end
 end
\ No newline at end of file
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
index 833f17b..e97b908 100644
--- a/@BearImp/calcCap.m
+++ b/@BearImp/calcCap.m
@@ -15,7 +15,7 @@ assert(obj.up2date.T,'Schmierstoffparameter noch nicht berechnet')
 assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
 assert(obj.up2date.B,'Belastungsverteilung noch nicht berechnet')
 assert(obj.up2date.H,'Schmierfilmdicke noch nicht berechnet')
-L = obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; H = obj.H; AddOn = obj.AddOn; psi = obj.psi;
+L = obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; H = obj.H; AddOn = obj.AddOn; psi = obj.psi; T_Oil = obj.T_Oil;
 C = obj.C;
 C.method = obj.method.C;
 if ~isfield(C,'k_C')
@@ -69,13 +69,20 @@ for posBall_conductive=1:L.numberOfConductiveBalls
         
     contactInd=B.contactInd(posBall_conductive,:);
     
-    switch C.method.rhoRatio % all from Brüser 1972
-        case 'Hartung' % (4.33, p. 43) (default)
+    switch C.method.rhoRatio
+        case 'Hartung' % Brüser 1972 (4.33, p. 43) (default)
             C.rhoRatio = @(p) 1  +  (9.73e-3 .* (p.*0.101972e-6).^0.75)  ./  (T.nu_38.*1e6).^0.0385;
-        case 'Dowson'  % (4.27, p. 41)
+        case 'Dowson'  % Brüser 1972 (4.27, p. 41)
             C.rhoRatio = @(p) 1 + 0.0057.*(p.*0.101972e-6)./(1+0.0165.*(p.*0.101972e-6)); % assuming [p] = kP/mm²
-        case 'Bradbury'% (4.30, p. 43)
+        case 'Bradbury'% Brüser 1972 (4.30, p. 43)
             C.rhoRatio = @(p) 1 ./ (1 - 0.087.*log(0.0684.*(p.*0.101972e-6)+1));          % assuming [p] = kP/mm²
+        case 'Tait'    % Bair 2007 (Table 6.4, p. 123)
+            K_0prime = 11; K_00 = 9e9; beta_K = 6.5e-3; a_V = 8e-4;
+            K_0 = K_00 * exp(-beta_K*(T_Oil+273.15));    % by Fakhreddine and Zoller in Bair 2007 (4.12, p. 69)
+            VbyV_0 = @(p) 1 - 1/(1+K_0prime)*log(1+p./K_0*(1+K_0prime)); % Bair 2007 (4.6 , p. 68)
+            V_0byV_R = 1./(1 - S.alpha_rho.*(T_Oil-15));                 % Bair 2007 (4.13, p. 70)
+            %V_0byV_R = 1 + a_V*(T_Oil-15);                              % Bair 2007 (4.14, p.70) depreciated
+            C.rhoRatio = @(p) 1./(VbyV_0(p) .* V_0byV_R);
     end
     C.epsilon_primeOil = @(p) (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*C.rhoRatio(p))  ./  (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*C.rhoRatio(p));
     if C.method.pressureDistribution == 'on'
@@ -90,12 +97,18 @@ for posBall_conductive=1:L.numberOfConductiveBalls
     else
         C.C_Hertz(:,contactInd,posBall_conductive) = obj.epsilon_0 .* C.epsilon_primeOil(temp_p(:,contactInd)) .* temp_A(:,contactInd) ./ temp_h_0(:,contactInd);
     end
+
+    if C.method.roughnessCorrection == 'Napel'
+        A = (L.Rz_RE + [L.Rz_IR;L.Rz_OR])/4;
+        C.napelFactor = 1./sqrt(1 - (A./temp_h_0(:,contactInd)).^2);
+        C.C_Hertz(:,contactInd,posBall_conductive) = C.C_Hertz(:,contactInd,posBall_conductive) .* C.napelFactor;
+    end
  
     R_el_single(:,contactInd,posBall_conductive) = 1./C.k_R .* S.rho_el .* temp_h_0(:,contactInd) ./ temp_A(:,contactInd);
     C.DeltaR_i(contactInd)    = tempR_li - tempR_RE - temp_h_0(1,contactInd);
     C.DeltaR_a(contactInd)    = tempR_la - tempR_RE - temp_h_0(2,contactInd);
 
-    if size(B.noContactInd,2) % Wenn Kraftfreie Wälzkörper vorhanden
+    if size(B.noContactInd,2) % Wenn kraftfreie Wälzkörper vorhanden
         C.DeltaR_i(B.noContactInd(posBall_conductive,:)') = tempR_li - tempR_RE - temp_s(1,B.noContactInd(posBall_conductive,:)');
         C.DeltaR_a(B.noContactInd(posBall_conductive,:)') = tempR_la - tempR_RE - temp_s(2,B.noContactInd(posBall_conductive,:)');
     end
diff --git a/@BearImp/calcClear.m b/@BearImp/calcClear.m
index 1a32a75..2423ccf 100644
--- a/@BearImp/calcClear.m
+++ b/@BearImp/calcClear.m
@@ -7,6 +7,11 @@ function calcClear(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler, Julius van der Kuip
 
+if ~isscalar(obj)
+    executeAllObj(obj,@calcClear)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.L,'Lager nicht gesetzt')
 L = obj.L; F_r = obj.F_r;
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index 4a97acc..eeb45ce 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -12,6 +12,12 @@ arguments
     options.Delta_Tia = 5 % Temperaturerhöhung des Innenrings gegenüber des Außenrings
     options.alpha = []
 end
+
+if ~isscalar(obj)
+    executeAllObj(obj,@calcGeo)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.L,'Lager nicht gesetzt')
 assert(obj.up2date.R,'Lagerspiel noch nicht berechnet')
diff --git a/@BearImp/calcLub.m b/@BearImp/calcLub.m
index 27d1eab..29513a4 100644
--- a/@BearImp/calcLub.m
+++ b/@BearImp/calcLub.m
@@ -5,6 +5,11 @@ function calcLub(obj)
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler
 
+if ~isscalar(obj)
+    executeAllObj(obj,@calcLub)
+    return
+end
+
 %% Parameter prüfen
 assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
 S = obj.S; T_Oil = obj.T_Oil;
@@ -23,18 +28,26 @@ switch T.method
         T.eta_00 = T.eta_040 - T.alpha_eta.*40;
         T.alpha_nu = (S.nu_100 - S.nu_40)/60;
         T.nu_0 = S.nu_40 - T.alpha_nu.*40;
-        T.eta = @(T) T.eta_00 + T.alpha_eta.*T;
-        T.eta_0 = T.eta(T_Oil);
-        T.nu_38 = T.nu_0 + T.alpha_nu.*38;
+        T.eta = @(theta) T.eta_00 + T.alpha_eta.*theta;
+        T.nu = @(theta) T.nu_0 + T.alpha_nu.*theta;
+    case 'DIN 51563'
+        assert(~any(isnan([S.nu_40 S.nu_100])),'Viscosity values not given for the selected lubricant. Choose method.T = ''Vogel''')
+        W_40  = log10(log10(S.nu_40 *1e6 + 0.8));                         % DIN 51563 (1)
+        W_100 = log10(log10(S.nu_100*1e6 + 0.8));                         % DIN 51563 (2)
+        m = (W_40-W_100)/(log10(100+273.15)-log10(40+273.15));            % DIN 51563 (3)
+        W_x = @(theta) m*(log10(40+273.15) - log10(theta+273.15)) + W_40; % DIN 51563 (4)
+        T.nu = @(theta) (10.^(10.^W_x(theta)) - 0.8)*1e-6;
+        T.eta = @(theta) T.rho(theta).*nu(theta);
     case 'Vogel'
         assert(~any(isnan([S.B S.C S.K])),'Vogel-Parameters not given for the selected lubricant. Choose method.T = ''linear''')
-        T.eta = @(T) S.K .* exp(S.B ./ (T + S.C));
-        T.eta_0    = T.eta(T_Oil);
+        T.eta = @(theta) S.K .* exp(S.B ./ (theta + S.C));
         T.eta_040  = T.eta(  40 );
         T.eta_0100 = T.eta( 100 );
-        T.nu_38 = T.eta(38) ./ T.rho(38);
+        T.nu = @(theta) T.eta(theta) ./ T.rho(theta);
 end
 T.alpha_etaT = log(T.eta_040/T.eta_0100)/60;
+T.eta_0    = T.eta(T_Oil);
+T.nu_38 = T.nu(38);
 
 %% Attribute ändern
 obj.T = T;
diff --git a/@BearImp/set.m b/@BearImp/set.m
index e512d21..26c5a06 100644
--- a/@BearImp/set.m
+++ b/@BearImp/set.m
@@ -1,6 +1,6 @@
 function set(obj, name, value)
-%Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
-%auch mit BearImp-Arrays bzw. -Matrizen
+% Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
+% auch mit BearImp-Arrays bzw. -Matrizen
 
 %    pmd Berechnungstool Lagerimpedanz
 %    Autor: Steffen Puchtler
@@ -24,13 +24,23 @@ function set(obj, name, value)
 end
 
 function setObjArrayPropToSingleValue(objArray, name, value)
+    propStack = strsplit(name,'.');
     for ii = 1:numel(objArray)
-        objArray(ii).(name) = value;
+        switch numel(propStack)
+            case 1
+                objArray(ii).(name) = value;
+            case 2
+                objArray(ii).(propStack{1}).(propStack{2}) = value;
+            case 3
+                objArray(ii).(propStack{1}).(propStack{2}).(propStack{3}) = value;
+            otherwise
+                error('Deeper property stacks are not yet implemented')
+        end
     end
 end
 
 function setObjArrayPropToValueArray(objArray, name, valueArray)
     for ii = 1:numel(objArray)
-        objArray(ii).(name) = valueArray(ii);
+        setObjArrayPropToSingleValue(objArray(ii), name, valueArray(ii))
     end
 end
\ No newline at end of file
diff --git a/BearImpOptions.m b/BearImpOptions.m
index cb28646..5c3f2f7 100644
--- a/BearImpOptions.m
+++ b/BearImpOptions.m
@@ -29,7 +29,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.T = {'linear','Vogel','DIN 51563'};
             possibleOption.B = {'static','dynamic'};
             possibleOption.H.filmThicknessFormula = {'Hamrock/Dowson','Moes'};
             possibleOption.H.thermalCorrection = {'on','off'};
@@ -39,7 +39,8 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             possibleOption.C.outsideArea = {'k-factor','neglect','stateOfTheArt','Schneider_k_c','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D','Puchtler2025'};
             possibleOption.C.k_vh_factor = {'on','off'};
             possibleOption.C.pressureDistribution = {'on','off'};
-            possibleOption.C.rhoRatio = {'Hartung','Dowson','Bradbury'};
+            possibleOption.C.roughnessCorrection = {'off','Napel'};
+            possibleOption.C.rhoRatio = {'Hartung','Dowson','Bradbury','Tait'};
         end
 
         function defaultOption = Default
@@ -55,6 +56,7 @@ classdef BearImpOptions < handle & dynamicprops & matlab.mixin.CustomDisplay & m
             defaultOption.C.outsideArea = 'semianalytisch3D';
             defaultOption.C.k_vh_factor = 'off';
             defaultOption.C.pressureDistribution = 'on';
+            defaultOption.C.roughnessCorrection = 'off';
             defaultOption.C.rhoRatio = 'Hartung';
         end
     end
-- 
GitLab