diff --git a/@BearImp/BearImp.m b/@BearImp/BearImp.m
index dbf7c35900faa15ff9962f0a19cd163bfc5cdc7f..9e51a353cb7eb2da5ce48e480e807908bd0a232c 100644
--- a/@BearImp/BearImp.m
+++ b/@BearImp/BearImp.m
@@ -1,411 +1,411 @@
-classdef BearImp < handle & matlab.mixin.Copyable
-% pmd Berechnungstool Lagerimpedanz
-% Fachgebiet für Produktentwicklung und Maschinenelemente, TU Darmstadt
-% Autoren: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
-
- properties (Dependent, Access = public)
- % Eingangsparameter Zugriff
- F_r (1,1) double {mustBeNonnegative} % Radialkraft in N
- F_a (1,1) double {mustBeNonnegative} % Axialkraft in N
- omega (1,1) double {mustBePositive} % Winkelgeschwindigkeit in 1/s
- T_Oil (1,1) double {mustBeReal} % Öltemperatur in °C
- psi (1,:) double {mustBeReal} % Winkel, die vom Benutzer für die Kapazitätsberechnung verwendet wird, in rad
- resolutionPerRotation (1,1) double {mustBeNonnegative, mustBeInteger} % Anzahl an Zwischenschritten des Winkels psi innerhalb einer 360°-Drehung
- L (1,1) struct % Lagerparameter (gegeben)
- S (1,1) struct % Schmierstoffparameter (gegeben)
- method (1,1) struct % struct der Rechenmethoden mit möglichen Feldern T,R,G,B,H,Z
- AddOn (1,1) struct % Installierte AddOns
-
-
- % Structs mit Berechnungsparamtern
- T (1,1) struct % 2.1 Tribologie (berechnete Schmierstoffparameter) f( S T_Oil )
- R (1,1) struct % 2.2 Radiallagerluft f(L F_r )
- G (1,1) struct % 2.3 Geometrie (berechnet in Abhängigkeit des Lagerspiels) f(L, R T_Oil )
- B (1,1) struct % 2.4 Belastungsverteilung f(L, G F_r,psi,AddOn psi_calc psi_all)
- H (1,1) struct % 2.5 Schmierfilmdicke f(L,S,T, G,B T_Oil, omega )
- C (1,1) struct % 2.6 Kapazität und Widerstand f(L,S,T, G,B,H psi,AddOn psi_calc psi_all)
- Z (1,1) struct % 2.7 Impedanz f(L C )
- end
-
-
- properties (Access = private)
- % Eingangsparameter Speicher
- privateInputParameters = struct('F_r',nan,'F_a',nan,'omega',nan,'T_Oil',nan,'psi',nan,'psi_calc',nan,'psi_all',nan,'ind_psi_all',nan,'f',nan,'L',struct,'S',struct,'resolutionPerRotation',nan)
- privateMethod = struct
- privateResults = struct('T',struct,'R',struct,'G',struct,'B',struct,'H',struct,'C',struct,'Z',struct)
- privateAddOn = struct('Parallel_Computing_Toolbox', 'auto', 'Optimization_Toolbox', 'auto')
- end
-
- properties (Access = public)
- % gibt an, ob die Werte vorhanden und aktuell sind
- up2date = struct('L',false,'S',false,'T',false,'R',false,'G',false,'B',false,'H',false,'C',false,'Z',false)
- materials
- inputDataTableFile = 'InputData.xlsx';
- PlotLabelTableFile = 'PlotLabel.csv';
- end
-
- properties (Dependent, SetAccess = protected, GetAccess = public)
- psi_calc (1,:) double {mustBeReal} % Winkel, die benötigt werden, um die Kapazität berechnen zu können, (unter ausnutzung der Symmetrie und ohne Dopplungen) in rad
- psi_all (1,:) double {mustBeReal} % alle Winkel, die unter Ausnutzung der Symmetrie mit psi_calc berechnet wurden
- ind_psi_all (:,:) double {mustBeNonnegative, mustBeInteger} % Indizes überführen psi_calc zu psi_all (damit kann zB. h_min(ind_psi_all) zu einem maximalen Vektor formiert werden)
- end
-
- properties (Dependent, SetAccess = immutable, GetAccess = public)
- allUp2date (1,1) logical % true, wenn alle berechneten Größen auf dem aktuellen Stand sind
- ready2calc (1,1) logical % true, wenn alle Prameter gegeben zur Berechnung
- end
-
- properties (Constant)
- epsilon_0 = 8.854187817e-12; % Permittivität des Vakuums
- end
-
- methods (Access = public)
- function obj = BearImp(setup)
- % Konstruktor tmp ###########
- % BearImp -> erzeugt leeres Objekt
- % BearImp('default') -> erzeugt Objekt mit default Werten (s.u.)
- if nargin > 0
- switch setup
- case 'default'
- obj.setBearing('6205-C3');
- obj.setLube('FVA Referenz-Öl III');
- obj.F_r = 1000;
- obj.F_a = 0;
- obj.omega = 3000/60*2*pi;
- obj.T_Oil = 50;
- obj.resolutionPerRotation = 360;
- end
- end
-
- end
- end
-
- methods (Access = public) % Später protected? #######################
- % Berechnungsmethoden in einzelnen Dateien
- calcLub (obj)
- calcClear(obj)
- calcGeo (obj)
- calcLoad (obj)
- varargout = calcFilm(obj, options)
- calcCap (obj)
- calcImp (obj)
- end
- methods (Access = public)
- setBearing(obj,name)
- setFitting(obj,shaft,housing)
- setLube (obj,name)
- calculate (obj)
- calculateContinuous(obj,n)
- set(obj,name,value)
- value = get(obj,name)
- obj = copyToArray(obj,varargin)
- end
-
- methods
- % set-Methoden
- function set.F_r(obj,F_r)
- obj.up2date.B = false;
- obj.up2date.R = false;
- obj.privateInputParameters.F_r = F_r;
- end
- function set.F_a(obj,F_a)
- obj.up2date.B = false;
- obj.up2date.G = false;
- obj.privateInputParameters.F_a = F_a;
- end
- function set.omega(obj,omega)
- obj.up2date.H = false;
- obj.privateInputParameters.omega = omega;
- end
- function set.T_Oil(obj,T_Oil)
- obj.up2date.H = false;
- obj.up2date.T = false;
- obj.up2date.G = false;
- obj.privateInputParameters.T_Oil = T_Oil;
- end
- function set.resolutionPerRotation(obj,resolutionPerRotation)
- obj.up2date.B = false;
- obj.up2date.C = false;
- obj.privateInputParameters.resolutionPerRotation = resolutionPerRotation;
- obj.privateInputParameters.psi = (0:(2*pi)/(resolutionPerRotation):pi*2*(resolutionPerRotation-1)/resolutionPerRotation);
- end
- function set.psi(obj,psi)
- obj.up2date.C = false;
- obj.up2date.B = false;
- obj.privateInputParameters.psi = psi;
- obj.privateInputParameters.resolutionPerRotation = nan;
- end
- function set.psi_calc(obj,psi_calc)
- obj.up2date.C = false;
- obj.up2date.B = false;
- obj.privateInputParameters.psi_calc = psi_calc;
- obj.privateInputParameters.resolutionPerRotation = nan;
- end
- function set.psi_all(obj,psi_all)
- obj.up2date.C = false;
- obj.up2date.B = false;
- obj.privateInputParameters.psi_all = psi_all;
- obj.privateInputParameters.resolutionPerRotation = nan;
- end
- function set.ind_psi_all(obj,ind_psi_all)
-% obj.up2date.R = false;
-% obj.up2date.B = false;
- obj.privateInputParameters.ind_psi_all = ind_psi_all;
- end
- function set.AddOn(obj,AddOn)
- obj.up2date.B = false;
- obj.up2date.C = false;
-
- pct = AddOn.Parallel_Computing_Toolbox;
- switch pct
- case true
- pct = 'on';
- case false
- pct = 'off';
- otherwise
- end
-
- ot = AddOn.Optimization_Toolbox;
- switch ot
- case true
- ot = 'on';
- case false
- ot = 'off';
- otherwise
- end
-
- check = {'auto', 'on', 'off'};
- if ~ischar(ot)||~ischar(pct)
- error("Wrong Input for AddOn. You have to decide between 'on', 'off' and 'auto'.")
- elseif ~any(strcmp(pct,check)==1) || ~any(strcmp(ot,check)== 1)
- error("Check if you spelled the AddOn Input the right way! ('on' 'off' 'auto')")
- end
-
- AddOn.Optimization_Toolbox = ot;
- AddOn.Parallel_Computing_Toolbox = pct;
- obj.privateAddOn = AddOn;
- end
- function set.L(obj,L)
- obj.up2date.L = true;
- obj.up2date.R = false;
- obj.up2date.G = false;
- obj.up2date.B = false;
- obj.up2date.Z = false;
- obj.privateInputParameters.L = L;
- end
- function set.S(obj,S)
- obj.up2date.S = true;
- obj.up2date.T = false;
- obj.up2date.H = false;
- obj.up2date.C = false;
- obj.privateInputParameters.S = S;
- end
- function set.method(obj,method)
- possibleMethods.check(method,true);
- somethingChanged = false;
- for ii = {'T','R','G','B','H','C','Z'}
- if isfield(obj.method,ii) && isfield (method,ii)
- if ~strcmp(obj.method.(ii{1}), method.(ii{1}))
- obj.up2date.(ii{1}) = false;
- somethingChanged = true;
- end
- elseif isfield(obj.method,ii{1}) || isfield (method,ii{1})
- obj.up2date.(ii{1}) = false;
- somethingChanged = true;
- end
- end
- if somethingChanged
- obj.privateMethod = method;
- else
- warning('Keine Methode geändert')
- end
- end
-
- function set.up2date(obj,up2date)
- try % Auf integrität prüfen
- assert(isstruct(up2date))
- vals = struct2array(up2date);
- assert(all(islogical(vals)))
- catch ME
- cause = MException('MATLAB:Lagerimpedanz:up2dateLogicError','up2date darf nur Variablen vom Typ logical enthalten');
- ME = addCause(ME,cause);
- rethrow(ME)
- end
-
- % Abhängige Variablen auf nicht up2date setzen
- if ~up2date.T, up2date.H = false; up2date.C = false; end
- if ~up2date.R, up2date.G = false; end
- if ~up2date.G, up2date.B = false; up2date.H = false; up2date.C = false; end
- if ~up2date.B, up2date.H = false; up2date.C = false; end
- if ~up2date.H, up2date.C = false; end
- if ~up2date.C, up2date.Z = false; end
-
- obj.up2date = up2date;
- end
- function set.T(obj,T)
- obj.up2date.H = false;
- obj.up2date.C = false;
- obj.privateResults.T = T;
- end
- function set.R(obj,R)
- obj.up2date.G = false;
- obj.privateResults.R = R;
- end
- function set.G(obj,G)
- obj.up2date.B = false;
- obj.up2date.H = false;
- obj.up2date.C = false;
- obj.privateResults.G = G;
- end
- function set.B(obj,B)
- obj.up2date.H = false;
- obj.up2date.C = false;
- obj.privateResults.B = B;
- end
- function set.H(obj,H)
- obj.up2date.C = false;
- obj.privateResults.H = H;
- end
- function set.C(obj,C)
- obj.up2date.Z = false;
- obj.privateResults.C = C;
- end
- function set.Z(obj,Z)
- obj.privateResults.Z = Z;
- end
- end
-
- methods
- % get-Methoden
- function F_r = get.F_r(obj)
- F_r = obj.privateInputParameters.F_r;
- end
- function F_a = get.F_a(obj)
- F_a = obj.privateInputParameters.F_a;
- end
- function omega = get.omega(obj)
- omega = obj.privateInputParameters.omega;
- end
- function T_Oil = get.T_Oil(obj)
- T_Oil = obj.privateInputParameters.T_Oil;
- end
- function psi = get.psi(obj)
- psi = obj.privateInputParameters.psi;
- end
- function psi_calc = get.psi_calc(obj)
- psi_calc = obj.privateInputParameters.psi_calc;
- end
- function psi_all = get.psi_all(obj)
- psi_all = obj.privateInputParameters.psi_all;
- end
- function ind_psi_all = get.ind_psi_all(obj)
- ind_psi_all = obj.privateInputParameters.ind_psi_all;
- end
- function resolutionPerRotation = get.resolutionPerRotation(obj)
- resolutionPerRotation = obj.privateInputParameters.resolutionPerRotation;
- end
- function AddOn = get.AddOn(obj)
-
- pct = obj.privateAddOn.Parallel_Computing_Toolbox;
- ot = obj.privateAddOn.Optimization_Toolbox;
-
- switch pct
- case 'on'
- AddOn.Parallel_Computing_Toolbox = true;
- case 'off'
- AddOn.Parallel_Computing_Toolbox = false;
- case 'auto'
- table_AddOns = matlab.addons.installedAddons;
- table_AddOns = table2array(table_AddOns(:,1));
-
- if any(strcmp(table_AddOns, "Parallel Computing Toolbox"),1) && matlab.addons.isAddonEnabled("Parallel Computing Toolbox")
- AddOn.Parallel_Computing_Toolbox = true;
- else
- AddOn.Parallel_Computing_Toolbox = false;
- disp("To speed up the calculation you can install the Parallel_Computing_Toolbox")
- end
-
- end
-
- switch ot
- case 'on'
- AddOn.Optimization_Toolbox = true;
- case 'off'
- AddOn.Optimization_Toolbox = false;
- case 'auto'
- table_AddOns = matlab.addons.installedAddons;
- table_AddOns = table2array(table_AddOns(:,1));
-
- if any(strcmp(table_AddOns, "Optimization Toolbox"),1) && matlab.addons.isAddonEnabled("Optimization Toolbox")
- AddOn.Optimization_Toolbox = true;
- else
- AddOn.Optimization_Toolbox = false;
- disp("To get more precise and faster calculations you can install the Optimization_Toolbox")
- end
-
- end
-
-
- end
- function L = get.L(obj)
- L = obj.privateInputParameters.L;
- end
- function S = get.S(obj)
- S = obj.privateInputParameters.S;
- end
- function method = get.method(obj)
- method = obj.privateMethod;
- end
- function T = get.T(obj)
- T = obj.privateResults.T;
- end
- function R = get.R(obj)
- R = obj.privateResults.R;
- end
- function G = get.G(obj)
- G = obj.privateResults.G;
- end
- function B = get.B(obj)
- B = obj.privateResults.B;
- end
- function H = get.H(obj)
- H = obj.privateResults.H;
- end
- function C = get.C(obj)
- C = obj.privateResults.C;
- end
- function Z = get.Z(obj)
- Z = obj.privateResults.Z;
- end
- function allUp2date = get.allUp2date(obj)
- allUp2date = all(struct2array(obj.up2date));
- end
- function ready2calc = get.ready2calc(obj)
- ready2calc = all(~isnan([obj.F_r obj.F_a obj.omega obj.T_Oil obj.psi])) &&...
- obj.up2date.L &&...
- obj.up2date.S;
- end
- function materials = get.materials(obj)
- if ~isstruct(obj.materials)
- materialTable = readtable(obj.inputDataTableFile,'Sheet','Materials');
- materialTable = removevars(materialTable,'possible_el_behaviours');
- for ii = 1:height(materialTable)
- materials.(materialTable(ii,:).Material{:}) = table2struct(materialTable(ii,:));
- end
- obj.materials = materials;
- else
- materials = obj.materials;
- end
- end
- end
-
- methods (Access = protected)
- function executeAllObj(obj,functionHandle)
- % Führt eine Funktion für alle Einträge von obj aus
- % Funktioniert aktuell nur für Funktionen ohne Übergabeparameter
- for ii = 1:numel(obj)
- feval(functionHandle,obj(ii))
- end
- end
- end
+classdef BearImp < handle & matlab.mixin.Copyable
+% pmd Berechnungstool Lagerimpedanz
+% Fachgebiet für Produktentwicklung und Maschinenelemente, TU Darmstadt
+% Autoren: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
+
+ properties (Dependent, Access = public)
+ % Eingangsparameter Zugriff
+ F_r (1,1) double {mustBeNonnegative} % Radialkraft in N
+ F_a (1,1) double {mustBeNonnegative} % Axialkraft in N
+ omega (1,1) double {mustBePositive} % Winkelgeschwindigkeit in 1/s
+ T_Oil (1,1) double {mustBeReal} % Öltemperatur in °C
+ psi (1,:) double {mustBeReal} % Winkel, die vom Benutzer für die Kapazitätsberechnung verwendet wird, in rad
+ resolutionPerRotation (1,1) double {mustBeNonnegative, mustBeInteger} % Anzahl an Zwischenschritten des Winkels psi innerhalb einer 360°-Drehung
+ L (1,1) struct % Lagerparameter (gegeben)
+ S (1,1) struct % Schmierstoffparameter (gegeben)
+ method (1,1) struct % struct der Rechenmethoden mit möglichen Feldern T,R,G,B,H,Z
+ AddOn (1,1) struct % Installierte AddOns
+
+
+ % Structs mit Berechnungsparamtern
+ T (1,1) struct % 2.1 Tribologie (berechnete Schmierstoffparameter) f( S T_Oil )
+ R (1,1) struct % 2.2 Radiallagerluft f(L F_r )
+ G (1,1) struct % 2.3 Geometrie (berechnet in Abhängigkeit des Lagerspiels) f(L, R T_Oil )
+ B (1,1) struct % 2.4 Belastungsverteilung f(L, G F_r,psi,AddOn psi_calc psi_all)
+ H (1,1) struct % 2.5 Schmierfilmdicke f(L,S,T, G,B T_Oil, omega )
+ C (1,1) struct % 2.6 Kapazität und Widerstand f(L,S,T, G,B,H psi,AddOn psi_calc psi_all)
+ Z (1,1) struct % 2.7 Impedanz f(L C )
+ end
+
+
+ properties (Access = private)
+ % Eingangsparameter Speicher
+ privateInputParameters = struct('F_r',nan,'F_a',nan,'omega',nan,'T_Oil',nan,'psi',nan,'psi_calc',nan,'psi_all',nan,'ind_psi_all',nan,'f',nan,'L',struct,'S',struct,'resolutionPerRotation',nan)
+ privateMethod = struct
+ privateResults = struct('T',struct,'R',struct,'G',struct,'B',struct,'H',struct,'C',struct,'Z',struct)
+ privateAddOn = struct('Parallel_Computing_Toolbox', 'auto', 'Optimization_Toolbox', 'auto')
+ end
+
+ properties (Access = public)
+ % gibt an, ob die Werte vorhanden und aktuell sind
+ up2date = struct('L',false,'S',false,'T',false,'R',false,'G',false,'B',false,'H',false,'C',false,'Z',false)
+ materials
+ inputDataTableFile = 'InputData.xlsx';
+ PlotLabelTableFile = 'PlotLabel.csv';
+ end
+
+ properties (Dependent, SetAccess = protected, GetAccess = public)
+ psi_calc (1,:) double {mustBeReal} % Winkel, die benötigt werden, um die Kapazität berechnen zu können, (unter ausnutzung der Symmetrie und ohne Dopplungen) in rad
+ psi_all (1,:) double {mustBeReal} % alle Winkel, die unter Ausnutzung der Symmetrie mit psi_calc berechnet wurden
+ ind_psi_all (:,:) double {mustBeNonnegative, mustBeInteger} % Indizes überführen psi_calc zu psi_all (damit kann zB. h_min(ind_psi_all) zu einem maximalen Vektor formiert werden)
+ end
+
+ properties (Dependent, SetAccess = immutable, GetAccess = public)
+ allUp2date (1,1) logical % true, wenn alle berechneten Größen auf dem aktuellen Stand sind
+ ready2calc (1,1) logical % true, wenn alle Prameter gegeben zur Berechnung
+ end
+
+ properties (Constant)
+ epsilon_0 = 8.854187817e-12; % Permittivität des Vakuums
+ end
+
+ methods (Access = public)
+ function obj = BearImp(setup)
+ % Konstruktor tmp ###########
+ % BearImp -> erzeugt leeres Objekt
+ % BearImp('default') -> erzeugt Objekt mit default Werten (s.u.)
+ if nargin > 0
+ switch setup
+ case 'default'
+ obj.setBearing('6205-C3');
+ obj.setLube('FVA Referenz-Öl III');
+ obj.F_r = 1000;
+ obj.F_a = 0;
+ obj.omega = 3000/60*2*pi;
+ obj.T_Oil = 50;
+ obj.resolutionPerRotation = 360;
+ end
+ end
+
+ end
+ end
+
+ methods (Access = public) % Später protected? #######################
+ % Berechnungsmethoden in einzelnen Dateien
+ calcLub (obj)
+ calcClear(obj)
+ calcGeo (obj)
+ calcLoad (obj)
+ varargout = calcFilm(obj, options)
+ calcCap (obj)
+ calcImp (obj)
+ end
+ methods (Access = public)
+ setBearing(obj,name)
+ setFitting(obj,shaft,housing)
+ setLube (obj,name)
+ calculate (obj)
+ calculateContinuous(obj,n)
+ set(obj,name,value)
+ value = get(obj,name)
+ obj = copyToArray(obj,varargin)
+ end
+
+ methods
+ % set-Methoden
+ function set.F_r(obj,F_r)
+ obj.up2date.B = false;
+ obj.up2date.R = false;
+ obj.privateInputParameters.F_r = F_r;
+ end
+ function set.F_a(obj,F_a)
+ obj.up2date.B = false;
+ obj.up2date.G = false;
+ obj.privateInputParameters.F_a = F_a;
+ end
+ function set.omega(obj,omega)
+ obj.up2date.H = false;
+ obj.privateInputParameters.omega = omega;
+ end
+ function set.T_Oil(obj,T_Oil)
+ obj.up2date.H = false;
+ obj.up2date.T = false;
+ obj.up2date.G = false;
+ obj.privateInputParameters.T_Oil = T_Oil;
+ end
+ function set.resolutionPerRotation(obj,resolutionPerRotation)
+ obj.up2date.B = false;
+ obj.up2date.C = false;
+ obj.privateInputParameters.resolutionPerRotation = resolutionPerRotation;
+ obj.privateInputParameters.psi = (0:(2*pi)/(resolutionPerRotation):pi*2*(resolutionPerRotation-1)/resolutionPerRotation);
+ end
+ function set.psi(obj,psi)
+ obj.up2date.C = false;
+ obj.up2date.B = false;
+ obj.privateInputParameters.psi = psi;
+ obj.privateInputParameters.resolutionPerRotation = nan;
+ end
+ function set.psi_calc(obj,psi_calc)
+ obj.up2date.C = false;
+ obj.up2date.B = false;
+ obj.privateInputParameters.psi_calc = psi_calc;
+ obj.privateInputParameters.resolutionPerRotation = nan;
+ end
+ function set.psi_all(obj,psi_all)
+ obj.up2date.C = false;
+ obj.up2date.B = false;
+ obj.privateInputParameters.psi_all = psi_all;
+ obj.privateInputParameters.resolutionPerRotation = nan;
+ end
+ function set.ind_psi_all(obj,ind_psi_all)
+% obj.up2date.R = false;
+% obj.up2date.B = false;
+ obj.privateInputParameters.ind_psi_all = ind_psi_all;
+ end
+ function set.AddOn(obj,AddOn)
+ obj.up2date.B = false;
+ obj.up2date.C = false;
+
+ pct = AddOn.Parallel_Computing_Toolbox;
+ switch pct
+ case true
+ pct = 'on';
+ case false
+ pct = 'off';
+ otherwise
+ end
+
+ ot = AddOn.Optimization_Toolbox;
+ switch ot
+ case true
+ ot = 'on';
+ case false
+ ot = 'off';
+ otherwise
+ end
+
+ check = {'auto', 'on', 'off'};
+ if ~ischar(ot)||~ischar(pct)
+ error("Wrong Input for AddOn. You have to decide between 'on', 'off' and 'auto'.")
+ elseif ~any(strcmp(pct,check)==1) || ~any(strcmp(ot,check)== 1)
+ error("Check if you spelled the AddOn Input the right way! ('on' 'off' 'auto')")
+ end
+
+ AddOn.Optimization_Toolbox = ot;
+ AddOn.Parallel_Computing_Toolbox = pct;
+ obj.privateAddOn = AddOn;
+ end
+ function set.L(obj,L)
+ obj.up2date.L = true;
+ obj.up2date.R = false;
+ obj.up2date.G = false;
+ obj.up2date.B = false;
+ obj.up2date.Z = false;
+ obj.privateInputParameters.L = L;
+ end
+ function set.S(obj,S)
+ obj.up2date.S = true;
+ obj.up2date.T = false;
+ obj.up2date.H = false;
+ obj.up2date.C = false;
+ obj.privateInputParameters.S = S;
+ end
+ function set.method(obj,method)
+ possibleMethods.check(method,true);
+ somethingChanged = false;
+ for ii = {'T','R','G','B','H','C','Z'}
+ if isfield(obj.method,ii) && isfield (method,ii)
+ if ~strcmp(obj.method.(ii{1}), method.(ii{1}))
+ obj.up2date.(ii{1}) = false;
+ somethingChanged = true;
+ end
+ elseif isfield(obj.method,ii{1}) || isfield (method,ii{1})
+ obj.up2date.(ii{1}) = false;
+ somethingChanged = true;
+ end
+ end
+ if somethingChanged
+ obj.privateMethod = method;
+ else
+ warning('Keine Methode geändert')
+ end
+ end
+
+ function set.up2date(obj,up2date)
+ try % Auf integrität prüfen
+ assert(isstruct(up2date))
+ vals = struct2array(up2date);
+ assert(all(islogical(vals)))
+ catch ME
+ cause = MException('MATLAB:Lagerimpedanz:up2dateLogicError','up2date darf nur Variablen vom Typ logical enthalten');
+ ME = addCause(ME,cause);
+ rethrow(ME)
+ end
+
+ % Abhängige Variablen auf nicht up2date setzen
+ if ~up2date.T, up2date.H = false; up2date.C = false; end
+ if ~up2date.R, up2date.G = false; end
+ if ~up2date.G, up2date.B = false; up2date.H = false; up2date.C = false; end
+ if ~up2date.B, up2date.H = false; up2date.C = false; end
+ if ~up2date.H, up2date.C = false; end
+ if ~up2date.C, up2date.Z = false; end
+
+ obj.up2date = up2date;
+ end
+ function set.T(obj,T)
+ obj.up2date.H = false;
+ obj.up2date.C = false;
+ obj.privateResults.T = T;
+ end
+ function set.R(obj,R)
+ obj.up2date.G = false;
+ obj.privateResults.R = R;
+ end
+ function set.G(obj,G)
+ obj.up2date.B = false;
+ obj.up2date.H = false;
+ obj.up2date.C = false;
+ obj.privateResults.G = G;
+ end
+ function set.B(obj,B)
+ obj.up2date.H = false;
+ obj.up2date.C = false;
+ obj.privateResults.B = B;
+ end
+ function set.H(obj,H)
+ obj.up2date.C = false;
+ obj.privateResults.H = H;
+ end
+ function set.C(obj,C)
+ obj.up2date.Z = false;
+ obj.privateResults.C = C;
+ end
+ function set.Z(obj,Z)
+ obj.privateResults.Z = Z;
+ end
+ end
+
+ methods
+ % get-Methoden
+ function F_r = get.F_r(obj)
+ F_r = obj.privateInputParameters.F_r;
+ end
+ function F_a = get.F_a(obj)
+ F_a = obj.privateInputParameters.F_a;
+ end
+ function omega = get.omega(obj)
+ omega = obj.privateInputParameters.omega;
+ end
+ function T_Oil = get.T_Oil(obj)
+ T_Oil = obj.privateInputParameters.T_Oil;
+ end
+ function psi = get.psi(obj)
+ psi = obj.privateInputParameters.psi;
+ end
+ function psi_calc = get.psi_calc(obj)
+ psi_calc = obj.privateInputParameters.psi_calc;
+ end
+ function psi_all = get.psi_all(obj)
+ psi_all = obj.privateInputParameters.psi_all;
+ end
+ function ind_psi_all = get.ind_psi_all(obj)
+ ind_psi_all = obj.privateInputParameters.ind_psi_all;
+ end
+ function resolutionPerRotation = get.resolutionPerRotation(obj)
+ resolutionPerRotation = obj.privateInputParameters.resolutionPerRotation;
+ end
+ function AddOn = get.AddOn(obj)
+
+ pct = obj.privateAddOn.Parallel_Computing_Toolbox;
+ ot = obj.privateAddOn.Optimization_Toolbox;
+
+ switch pct
+ case 'on'
+ AddOn.Parallel_Computing_Toolbox = true;
+ case 'off'
+ AddOn.Parallel_Computing_Toolbox = false;
+ case 'auto'
+ table_AddOns = matlab.addons.installedAddons;
+ table_AddOns = table2array(table_AddOns(:,1));
+
+ if any(strcmp(table_AddOns, "Parallel Computing Toolbox"),1) && matlab.addons.isAddonEnabled("Parallel Computing Toolbox")
+ AddOn.Parallel_Computing_Toolbox = true;
+ else
+ AddOn.Parallel_Computing_Toolbox = false;
+ disp("To speed up the calculation you can install the Parallel_Computing_Toolbox")
+ end
+
+ end
+
+ switch ot
+ case 'on'
+ AddOn.Optimization_Toolbox = true;
+ case 'off'
+ AddOn.Optimization_Toolbox = false;
+ case 'auto'
+ table_AddOns = matlab.addons.installedAddons;
+ table_AddOns = table2array(table_AddOns(:,1));
+
+ if any(strcmp(table_AddOns, "Optimization Toolbox"),1) && matlab.addons.isAddonEnabled("Optimization Toolbox")
+ AddOn.Optimization_Toolbox = true;
+ else
+ AddOn.Optimization_Toolbox = false;
+ disp("To get more precise and faster calculations you can install the Optimization_Toolbox")
+ end
+
+ end
+
+
+ end
+ function L = get.L(obj)
+ L = obj.privateInputParameters.L;
+ end
+ function S = get.S(obj)
+ S = obj.privateInputParameters.S;
+ end
+ function method = get.method(obj)
+ method = obj.privateMethod;
+ end
+ function T = get.T(obj)
+ T = obj.privateResults.T;
+ end
+ function R = get.R(obj)
+ R = obj.privateResults.R;
+ end
+ function G = get.G(obj)
+ G = obj.privateResults.G;
+ end
+ function B = get.B(obj)
+ B = obj.privateResults.B;
+ end
+ function H = get.H(obj)
+ H = obj.privateResults.H;
+ end
+ function C = get.C(obj)
+ C = obj.privateResults.C;
+ end
+ function Z = get.Z(obj)
+ Z = obj.privateResults.Z;
+ end
+ function allUp2date = get.allUp2date(obj)
+ allUp2date = all(struct2array(obj.up2date));
+ end
+ function ready2calc = get.ready2calc(obj)
+ ready2calc = all(~isnan([obj.F_r obj.F_a obj.omega obj.T_Oil obj.psi])) &&...
+ obj.up2date.L &&...
+ obj.up2date.S;
+ end
+ function materials = get.materials(obj)
+ if ~isstruct(obj.materials)
+ materialTable = readtable(obj.inputDataTableFile,'Sheet','Materials');
+ materialTable = removevars(materialTable,'possible_el_behaviours');
+ for ii = 1:height(materialTable)
+ materials.(materialTable(ii,:).Material{:}) = table2struct(materialTable(ii,:));
+ end
+ obj.materials = materials;
+ else
+ materials = obj.materials;
+ end
+ end
+ end
+
+ methods (Access = protected)
+ function executeAllObj(obj,functionHandle)
+ % Führt eine Funktion für alle Einträge von obj aus
+ % Funktioniert aktuell nur für Funktionen ohne Übergabeparameter
+ for ii = 1:numel(obj)
+ feval(functionHandle,obj(ii))
+ end
+ end
+ end
end
\ No newline at end of file
diff --git a/@BearImp/calcCap.m b/@BearImp/calcCap.m
index 05e2bac89842efec7ed8df893bf1a614f7f97046..672d6ee1a3ae13176701993b051769fd0222a58a 100644
--- a/@BearImp/calcCap.m
+++ b/@BearImp/calcCap.m
@@ -1,478 +1,478 @@
-function calcCap(obj)
-%2.6 Kapazität und Widerstand C = C(L,S,T,G,B,H)
-%Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
-%und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wird der
-%Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
-%berechnet.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
-
-%% Parameter prüfen
-assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
-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;
-C = struct;
-method = possibleMethods.addDefault(obj.method).C;
-if ~isfield(C,'k_C')
- C.k_C = 1;
-end
-if ~isfield(C,'k_R')
- C.k_R = 1;
-end
-
-%% Berechnung
-
-assert(~L.allBallsSameMaterial || L.RE_A.electric_behaviour == "conductor",'Die elektrischen Eigenschaften des Lagers können nur von RE berechnet werden, die ein leitendes Material besitzen. Bitte Material der RE prüfen.')
-
-if ~isnan(obj.psi_calc) % Falls psi_calc gegeben ist
- psi = obj.psi_calc; % psi_calc in calcCap als psi verwenden
-end
-
-B.numOfCagePositions = length(psi);
-
-if L.allBallsSameMaterial %|| strcmp(L.RE_order,'AAAAAAAAA')
- ballMaterialInd=1;
-else
- ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
-end
-
-C.C_Hertz = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
-C.phi_1i = nan (1,B.numOfCagePositions,L.numberOfConductiveBalls);
-C.phi_1a = nan (1,B.numOfCagePositions,L.numberOfConductiveBalls);
-C.C_out = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
-C.c_i = nan (1,B.numOfCagePositions);
-C.c_a = nan (1,B.numOfCagePositions);
-C.g_1i = zeros( B.numOfCagePositions,L.numberOfConductiveBalls);
-C.g_1a = zeros( B.numOfCagePositions,L.numberOfConductiveBalls);
-R_el_single = inf (2,B.numOfCagePositions,L.numberOfConductiveBalls);
-
-for posBall_conductive=1:L.numberOfConductiveBalls
-
- temp_s = B.s(:,:,posBall_conductive);
- temp_p = H.p(:,:,posBall_conductive);
- temp_A = H.A(:,:,posBall_conductive);
- temp_a = H.a(:,:,posBall_conductive);
- temp_b = H.b(:,:,posBall_conductive);
- temp_h_0 = H.h_0(:,:,posBall_conductive);
- tempR_li = G.R_li(ballMaterialInd(posBall_conductive));
- tempR_la = G.R_la(ballMaterialInd(posBall_conductive));
- tempR_RE = G.R_RE(ballMaterialInd(posBall_conductive));
-
- contactInd=B.contactInd(posBall_conductive,:);
-
- C.rhoRatio = 1 + (9.73e-3 .* (temp_p.*0.101972e-6).^0.75) ./ (T.nu_38.*1e6).^0.0385;
- C.epsilon_primeOil = (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*C.rhoRatio) ./ (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*C.rhoRatio);
- C.C_Hertz(:,contactInd,posBall_conductive) = obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ temp_h_0(:,contactInd);
-
- 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
- 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
-
- C.phi_1i(:,:,posBall_conductive) = temp_a(1,:)./G.D_RE(ballMaterialInd(posBall_conductive));
- C.phi_1a(:,:,posBall_conductive) = temp_a(2,:)./G.D_RE(ballMaterialInd(posBall_conductive));
- C.phi_2i = asin(L.B_i/G.D_RE(ballMaterialInd(posBall_conductive)));
- C.phi_2a = asin(L.B_a/G.D_RE(ballMaterialInd(posBall_conductive)));
-
- if any([strcmp(method.unloadedRE,{'Leander_Radial','LeanderSteffen'}) strcmp(method.outsideArea,{'Leander_Radial','LeanderSteffen'})])
- calcZg
- end
-
- if any(B.noContactInd,'all') % Wenn Kraftfreie Wälzkörper vorhanden
- switch method.unloadedRE
- case 'neglect'
- case 'Leander_Parallel'
- leander_parallel(B.noContactInd)
- case 'Leander_Radial'
- leander_radial(find(contactInd==0))
- case 'LeanderSteffen'
- leander_steffen(B.noContactInd)
- case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
- tobias_steffen(B.noContactInd)
- case 'semianalytisch3D'
-
- vecStruct=splitVec(contactInd,false);
-
- for recentVec=length(vecStruct)
- runSemianalytisch3D(vecStruct{recentVec})
- end
-
- end
- end
-
- switch method.outsideArea
- case 'neglect'
- case 'k-factor'
- C.C_Hertz(:,contactInd,posBall_conductive) = C.k_C .* obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ H.h_minth(:,contactInd);
- case 'Leander_Parallel'
- leander_parallel(B.contactInd)
- case 'Leander_Radial'
- leander_radial(find(contactInd==1))
- case 'LeanderSteffen'
- leander_steffen(B.contactInd)
- case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
- tobias_steffen(B.contactInd)
- case 'semianalytisch3D'
-
- vecStruct=splitVec(contactInd,true);
-
- for recentVec=1:length(vecStruct)
- runSemianalytisch3D(vecStruct{recentVec})
- end
- end
-
- C_el_single(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) + C.C_out(:,:,posBall_conductive);
-
-end
-
-if L.allBallsSameMaterial
-
- C_ind = check_convert_psi_calc2matrx(psi,L.Z,'useSymmetry');
-
- C.C_el_single = convert_vek2matr(C_el_single,C_ind);
- C.R_el_single = convert_vek2matr(R_el_single,C_ind);
-
- C_el_single = C.C_el_single;
- R_el_single = C.R_el_single;
-else
- check_convert_psi_calc2matrx(psi,L.Z,'noSymmetry');
- C.C_el_single = C_el_single;
- C.R_el_single = R_el_single;
-end
-
-C.C_el = 1 ./ (1./ sum(C_el_single(1,:,:),3) + 1./ sum( C_el_single(2,:,:) ,3));
-C.R_el = 1 ./ sum(1./ R_el_single(1,:,:),3) + 1./ sum(1./ R_el_single(2,:,:) ,3) ;
-
-%% Attribute ändern
-obj.C = C;
-obj.up2date.C = true;
-
-%% Helper-Functions
-function calcZg
-
- C.g_1i(contactInd,posBall_conductive) = asin( sin( C.phi_1i(1,contactInd,posBall_conductive)) .* (-C.DeltaR_i(:,contactInd)./tempR_li .* cos(C.phi_1i(1,contactInd,posBall_conductive)) + ...
- sqrt(1-C.DeltaR_i(:,contactInd).^2./tempR_li.^2 .* sin(C.phi_1i(1,contactInd,posBall_conductive)).^2)));
-
- C.g_1a(contactInd,posBall_conductive) = asin( sin( C.phi_1a(1,contactInd,posBall_conductive)) .* (-C.DeltaR_a(:,contactInd)./tempR_la .* cos(C.phi_1a(1,contactInd,posBall_conductive)) + ... %hier ein - statt +?
- sqrt(1-C.DeltaR_a(:,contactInd).^2./tempR_la.^2 .* sin(C.phi_1a(1,contactInd,posBall_conductive)).^2)));
-
- C.g_2i = acos(L.B_i/(2*tempR_li));
- C.g_2a = acos(L.B_a/(2*tempR_la));
-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);
-
- 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.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? ###
-
- 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.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;
- for ii = indices
- x1=double(C.g_1i(ii));
- x2=double(C.g_2i);
-
- C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_i(p)./C.h_i(p,t,ii), x1 , x2 , 0, pi/9);
- C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_a(p)./C.h_a(p,t,ii), C.g_1a(ii)+pi, C.g_2a+pi, 0, pi/9);
-% Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_i(p)./Z.h_i(p,t,ii), Z.g_1i(ii) , Z.g_2i , 0, pi/9);
-% Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_a(p)./Z.h_a(p,t,ii), Z.g_1a(ii)+pi, Z.g_2a+pi, 0, pi/9);
- end
-end
-
-function tobias_steffen(indices)
- C.h_i = @(alpha,beta,ii) (sqrt(tempR_li^2 - C.DeltaR_i(ii)^2*sin(alpha).^2) - C.DeltaR_i(ii)*cos(alpha)) ./ cos(beta) - tempR_RE;
- C.h_a = @(alpha,beta,ii) (sqrt(tempR_la^2 - C.DeltaR_a(ii)^2*sin(alpha).^2) - C.DeltaR_a(ii)*cos(alpha)) ./ cos(beta) - tempR_RE;
- if strcmp(method.unloadedRE,'TobiasSteffen_Kugelfläche')
- C.dA_i = @(~,beta,~) tempR_RE^2 * cos(beta);
- C.dA_a = C.dA_i;
- else
- C.dA_i = @(alpha,beta,ii) tempR_li./cos(beta) .* (1 - C.DeltaR_i(ii)*cos(alpha)./sqrt(tempR_li^2-C.DeltaR_i(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_i(alpha,beta,ii));
- C.dA_a = @(alpha,beta,ii) tempR_la./cos(beta) .* (1 - C.DeltaR_a(ii)*cos(alpha)./sqrt(tempR_la^2-C.DeltaR_a(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_a(alpha,beta,ii));
- end
- for ii = indices
- C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_i(alpha,beta,ii)./C.h_i(alpha,beta,ii), C.phi_1i(ii), C.phi_2i, 0, pi/3);
- C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_a(alpha,beta,ii)./C.h_a(alpha,beta,ii), C.phi_1a(ii), C.phi_2a, 0, pi/3);
- end
-end
-
-function runSemianalytisch3D(indices)
-
- if isfield(C,'C_out')
- C_semi_i = C.C_out(1,:,posBall_conductive);
- C_semi_a = C.C_out(2,:,posBall_conductive);
- else
- C_semi_i = nan( numel(temp_s)/2, 1);
- C_semi_a = nan( numel(temp_s)/2, 1);
- end
-
- 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);
-
- parfor runSemi = indices
- C_semi_i(runSemi) = feval(tmp_semi_i, runSemi); %da parfor loop nicht auf nested functions (semianalytisch3D) direkt zugreifen kann
- C_semi_a(runSemi )= feval(tmp_semi_a, runSemi);
- end
-
- 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);
- end
- end
-
- C.C_out(1,:,posBall_conductive)=C_semi_i;
- C.C_out(2,:,posBall_conductive)=C_semi_a;
-end
-
-function C_Zusatz=semianalytisch3D(s,R_WK,R_R,B_R,B,R_L,epsilon_r,varargin)
-% Semianalytische Berechnung der Kapazität von unbelasteten Wälzkörpern
-% sowie Randbereichen belasteter WK am Innen und Außenring
-% Autor: Steffen Puchtler
-% Herleitung: Masterthesis "Optimierung des Berechnungsverfahrens für die
-% elektrische Kapazität von EHD-Kontakten unter Berücksichtigung des realen
-% elektrischen Feldes", S. 23-26
-%
-% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r)
-% belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b)
-
-% Diskretisierung nicht Default:
-% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,0,0,nt,np)
-% belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b,nt,np)
-
-% Inputparameter:
-% S Minimaler Abstand zwischen Wälzkörper und Laufbahn m
-% R_WK Wälzkörperradius m
-% f Schmiegung (R_R / D_WK) -
-% B_R Rillenbreite m
-% B Breite des Lagers m
-% R_L Radius der Laufbahn m !!! R_L < 0 am Innenring
-% epsilon_r Relative Permittivität des Schmierstoffs -
-% a,b Halbachsen der Hertzschen Flächen m
-% nt,np Anzahl der Diskretisierngen in theta bzw. phi -
-
-% profile on
-
- switch nargin
- case 7
- nt = 60; np = 60;
- belastet = false;
- case 9
- nt = 110; np = 110;
- a = varargin{1};
- b = varargin{2};
- belastet = a~=0 && b~=0;
- case 11
- a = varargin{1};
- b = varargin{2};
- nt = varargin{3};
- np = varargin{4};
- belastet = a~=0 && b~=0;
- otherwise
- warning('The number of input arguments is wrong! Please check the input of semianalytisch3D.')
- end
-
- DeltaR = R_R-R_WK-s;
- R_T = R_L - R_R;
- Dz = R_R-sqrt(R_R^2-B_R^2/4);
- r_yz = @(phi) (R_WK-R_L+s).*cos(phi) + sqrt((R_L-Dz).^2-(R_WK-R_L+s).^2.*sin(phi).^2)*sign(R_L);
- if belastet
- Theta_0 = @(phi) a*sqrt(max(1/R_WK^2-phi.^2/b^2,0));
- else
- Theta_0 = @(~) 0;
- end
- Theta_1 = @(phi) atan(B_R/2./r_yz(phi));
- Theta_2 = @(phi) atan(B /2./r_yz(phi));
- if R_L > 0
- phi_1 = pi/2;
- else
- phi_1 = asin(R_L/(R_L-s-R_WK))*0.99;
- end
- p = linspace(0,phi_1,np);
- P = repmat(p',1,nt);
- T_t = zeros(size(P));
- for ii = 1:np
- T_t(ii,:) = linspace(Theta_0(p(ii)),Theta_1(p(ii)),nt);
- end
-
- r = zeros(size(T_t));
- for ii = 1:size(T_t,1)
- for jj = 1:size(T_t,2)
- theta = T_t(ii,jj);
- phi = P(ii,jj);
- fun = @(r) (R_T+R_R*sqrt(1-(r.*sin(theta)/R_R).^2)).^2 - (r.*sin(phi).*cos(theta)).^2 - (DeltaR+R_T+r.*cos(phi).*cos(theta)).^2;
-% options=optimset('Display','off','FunValCheck','off','OutputFcn',[],'PlotFcns',[],'TolX',eps);
- r(ii,jj) = fzero(fun,R_WK);
- end
- end
- clear phi theta
- invalid = r(:)>1;
- r2small = r(:)/R_WK<1.00125;
- if any(invalid)
- warning('%i Entries of r have been removed',sum(invalid))
- C_Zusatz=nan;
- return
- 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);
- firstInteg = zeros(1,np);
- for ii = 1:np
- firstInteg(ii) = trapz(T_t(ii,:),integrand_t(ii,:));
- end
- C_Rille = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
-
-
- T_k = zeros(size(P));
- for ii = 1:np
- T_k(ii,:) = linspace(Theta_1(p(ii)),Theta_2(p(ii)),nt);
- end
- integrand_k = R_WK^2./((r_yz(P)-s)./cos(T_k)-R_WK).*cos(T_k);
- firstInteg = zeros(1,np);
- for ii = 1:np
- firstInteg(ii) = trapz(T_k(ii,:),integrand_k(ii,:));
- end
- C_Rand = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
-
- C_Zusatz=C_Rille + C_Rand;
-% if nargout == 1
-% varargout = {C_Rille + C_Rand};
-% elseif nargout == 2
-% varargout = {C_Rille,C_Rand};
-% end
-
-% profile off
-% profile viewer
-
-end
-
-function seperateVecs=splitVec(vec,status)
- %da die parfor-Schleife nur Vektoren bearbeiten kann, die in Ganzzahl
- %Schritten auf- oder absteigen, wird in dieser Funktion jeder Vektor in
- %mehrere Vektoren aufgeteilt, die diese Bedingungen erfüllen.
-
- findInds = find(vec==status);
- startInds = unique([findInds(1) findInds(diff([0 findInds])>1)]); % Segment Start Indices
-
- if length(startInds)==1
- seperateVecs{1}=findInds;
- else
- findNonZero = [find(diff([0 findInds])>1)-1 length(findInds)];
- findNonZero(findNonZero == 0) = [];
- endInds = findInds(findNonZero); % Segment End Indices
-
- for runSection = 1:length(startInds)
- seperateVecs{runSection} = (startInds(runSection):endInds(runSection));
- end
- end
-end
-
-function ind_psi = check_convert_psi_calc2matrx(psi,Z,symCase)
- % Diese Funktion überprüft, ob der Input-Vektor psi alle Winkel
- % enthält, die benötigt werden, um die Gesamtkapazität des Wälzlagers
- % zu berechnen.
- % Außerdem wird eine Matrix ind_psi erstellt, die die benötigten
- % Indizes ausgibt, um dieberechneten Einzelkapazitäten in die richtige
- % Reihenfolge und Ordnung zu bringen.
-
- switch symCase
- case 'useSymmetry'
- symmetry = true;
- case 'noSymmetry'
- symmetry = false;
- otherwise
- error('Wrong entry for symCase! Possible Input: ''useSymmetry'' or ''noSymmetry''')
- end
- ind_psi = nan(obj.L.Z,length(obj.psi));
-
- psi_calc = psi;
-
- for run_C = 1 : length(obj.psi)
-
- psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
- psiNeedForIndCapa = round(psiNeedForIndCapa,8);
-
- while max(psiNeedForIndCapa) > (min(obj.psi) + 2*pi - 1e-8)
- psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) = psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
- end
-
- if ~isnan(obj.psi_calc)
-
- psi_check = psiNeedForIndCapa;
- symPoint = ceil(min(psi)/pi)*pi+pi;
- if symmetry %Es entsteht eine Symmetrie in der vertikalen Achse, wodurch die Berechnung von zB. der Kapazität nur für eine Hälfte der Umdrehung berechnet werden muss
- psi_check(psi_check > symPoint) = 2*symPoint - psi_check(psi_check > symPoint);
- tempPsiComp = 2*symPoint - psi_check(psi_check < symPoint - eps(1e8)); % um Rechenungenauigkeiten durch irrationale Zahl zu vermeiden, wird der Filten erweitert
-
- psi_check = uniquetol(psi_check, 1e-8);
- psi_check(ismembertol(psi_check,tempPsiComp,1e-8)) = []; %gespiegelte Punkte löschen
- end
- psi_raw = round(obj.psi,8);
- while (min(psi_check) < min(psi_raw)) || (max(psi_check) >= min(psi_raw) + 2*pi- eps(1e8))
- psi_check(psi_check<min(psi_raw)) = psi_check(psi_check<min(psi_raw))+ 2*pi;
- psi_check(psi_check>=min(psi_raw)+2*pi- eps(1e8)) = psi_check(psi_check>=min(psi_raw)+2*pi- eps(1e8))- 2*pi;
- end
-
- [~, ind_single_psi] = ismembertol(psi_check, psi_calc,1e-8);
- [~, ind_comp] = ismembertol(psiNeedForIndCapa, obj.psi_all,1e-8);
- ind_psi(:,run_C) = obj.ind_psi_all(ind_comp);
- else
- [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
- ind_psi(:,run_C) = ind_single_psi';
- end
- %Check, if all nessecary C_el_single are calculated
- assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!')
-
- end
-end
-
-function matrx_out = convert_vek2matr(A,B)
- % Diese Funktion erzeugt aus dem gewünschten Vektor und einer
- % Index-Matrix eine Matrix, gefüllt mit den Vektoreinträgen in der
- % Reihenfolge der Matrix-Indizes.
-
- matrx_out = nan(size(A,1),size(B,2),size(B,1));
- for run_1A = 1:size(A,1)
- for run_1B = 1:size(B,1)
- matrx_out(run_1A,:,run_1B) = A(run_1A,B(run_1B,:));
- end
- end
-end
+function calcCap(obj)
+%2.6 Kapazität und Widerstand C = C(L,S,T,G,B,H)
+%Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
+%und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wird der
+%Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
+%berechnet.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
+
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
+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;
+C = struct;
+method = possibleMethods.addDefault(obj.method).C;
+if ~isfield(C,'k_C')
+ C.k_C = 1;
+end
+if ~isfield(C,'k_R')
+ C.k_R = 1;
+end
+
+%% Berechnung
+
+assert(~L.allBallsSameMaterial || L.RE_A.electric_behaviour == "conductor",'Die elektrischen Eigenschaften des Lagers können nur von RE berechnet werden, die ein leitendes Material besitzen. Bitte Material der RE prüfen.')
+
+if ~isnan(obj.psi_calc) % Falls psi_calc gegeben ist
+ psi = obj.psi_calc; % psi_calc in calcCap als psi verwenden
+end
+
+B.numOfCagePositions = length(psi);
+
+if L.allBallsSameMaterial %|| strcmp(L.RE_order,'AAAAAAAAA')
+ ballMaterialInd=1;
+else
+ ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
+end
+
+C.C_Hertz = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.phi_1i = nan (1,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.phi_1a = nan (1,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.C_out = zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+C.c_i = nan (1,B.numOfCagePositions);
+C.c_a = nan (1,B.numOfCagePositions);
+C.g_1i = zeros( B.numOfCagePositions,L.numberOfConductiveBalls);
+C.g_1a = zeros( B.numOfCagePositions,L.numberOfConductiveBalls);
+R_el_single = inf (2,B.numOfCagePositions,L.numberOfConductiveBalls);
+
+for posBall_conductive=1:L.numberOfConductiveBalls
+
+ temp_s = B.s(:,:,posBall_conductive);
+ temp_p = H.p(:,:,posBall_conductive);
+ temp_A = H.A(:,:,posBall_conductive);
+ temp_a = H.a(:,:,posBall_conductive);
+ temp_b = H.b(:,:,posBall_conductive);
+ temp_h_0 = H.h_0(:,:,posBall_conductive);
+ tempR_li = G.R_li(ballMaterialInd(posBall_conductive));
+ tempR_la = G.R_la(ballMaterialInd(posBall_conductive));
+ tempR_RE = G.R_RE(ballMaterialInd(posBall_conductive));
+
+ contactInd=B.contactInd(posBall_conductive,:);
+
+ C.rhoRatio = 1 + (9.73e-3 .* (temp_p.*0.101972e-6).^0.75) ./ (T.nu_38.*1e6).^0.0385;
+ C.epsilon_primeOil = (S.epsilon_Oel + 2 + 2*(S.epsilon_Oel-1).*C.rhoRatio) ./ (S.epsilon_Oel + 2 - (S.epsilon_Oel-1).*C.rhoRatio);
+ C.C_Hertz(:,contactInd,posBall_conductive) = obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ temp_h_0(:,contactInd);
+
+ 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
+ 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
+
+ C.phi_1i(:,:,posBall_conductive) = temp_a(1,:)./G.D_RE(ballMaterialInd(posBall_conductive));
+ C.phi_1a(:,:,posBall_conductive) = temp_a(2,:)./G.D_RE(ballMaterialInd(posBall_conductive));
+ C.phi_2i = asin(L.B_i/G.D_RE(ballMaterialInd(posBall_conductive)));
+ C.phi_2a = asin(L.B_a/G.D_RE(ballMaterialInd(posBall_conductive)));
+
+ if any([strcmp(method.unloadedRE,{'Leander_Radial','LeanderSteffen'}) strcmp(method.outsideArea,{'Leander_Radial','LeanderSteffen'})])
+ calcZg
+ end
+
+ if any(B.noContactInd,'all') % Wenn Kraftfreie Wälzkörper vorhanden
+ switch method.unloadedRE
+ case 'neglect'
+ case 'Leander_Parallel'
+ leander_parallel(B.noContactInd)
+ case 'Leander_Radial'
+ leander_radial(find(contactInd==0))
+ case 'LeanderSteffen'
+ leander_steffen(B.noContactInd)
+ case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
+ tobias_steffen(B.noContactInd)
+ case 'semianalytisch3D'
+
+ vecStruct=splitVec(contactInd,false);
+
+ for recentVec=length(vecStruct)
+ runSemianalytisch3D(vecStruct{recentVec})
+ end
+
+ end
+ end
+
+ switch method.outsideArea
+ case 'neglect'
+ case 'k-factor'
+ C.C_Hertz(:,contactInd,posBall_conductive) = C.k_C .* obj.epsilon_0 .* C.epsilon_primeOil(:,contactInd) .* temp_A(:,contactInd) ./ H.h_minth(:,contactInd);
+ case 'Leander_Parallel'
+ leander_parallel(B.contactInd)
+ case 'Leander_Radial'
+ leander_radial(find(contactInd==1))
+ case 'LeanderSteffen'
+ leander_steffen(B.contactInd)
+ case {'TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche'}
+ tobias_steffen(B.contactInd)
+ case 'semianalytisch3D'
+
+ vecStruct=splitVec(contactInd,true);
+
+ for recentVec=1:length(vecStruct)
+ runSemianalytisch3D(vecStruct{recentVec})
+ end
+ end
+
+ C_el_single(:,:,posBall_conductive) = C.C_Hertz(:,:,posBall_conductive) + C.C_out(:,:,posBall_conductive);
+
+end
+
+if L.allBallsSameMaterial
+
+ C_ind = check_convert_psi_calc2matrx(psi,L.Z,'useSymmetry');
+
+ C.C_el_single = convert_vek2matr(C_el_single,C_ind);
+ C.R_el_single = convert_vek2matr(R_el_single,C_ind);
+
+ C_el_single = C.C_el_single;
+ R_el_single = C.R_el_single;
+else
+ check_convert_psi_calc2matrx(psi,L.Z,'noSymmetry');
+ C.C_el_single = C_el_single;
+ C.R_el_single = R_el_single;
+end
+
+C.C_el = 1 ./ (1./ sum(C_el_single(1,:,:),3) + 1./ sum( C_el_single(2,:,:) ,3));
+C.R_el = 1 ./ sum(1./ R_el_single(1,:,:),3) + 1./ sum(1./ R_el_single(2,:,:) ,3) ;
+
+%% Attribute ändern
+obj.C = C;
+obj.up2date.C = true;
+
+%% Helper-Functions
+function calcZg
+
+ C.g_1i(contactInd,posBall_conductive) = asin( sin( C.phi_1i(1,contactInd,posBall_conductive)) .* (-C.DeltaR_i(:,contactInd)./tempR_li .* cos(C.phi_1i(1,contactInd,posBall_conductive)) + ...
+ sqrt(1-C.DeltaR_i(:,contactInd).^2./tempR_li.^2 .* sin(C.phi_1i(1,contactInd,posBall_conductive)).^2)));
+
+ C.g_1a(contactInd,posBall_conductive) = asin( sin( C.phi_1a(1,contactInd,posBall_conductive)) .* (-C.DeltaR_a(:,contactInd)./tempR_la .* cos(C.phi_1a(1,contactInd,posBall_conductive)) + ... %hier ein - statt +?
+ sqrt(1-C.DeltaR_a(:,contactInd).^2./tempR_la.^2 .* sin(C.phi_1a(1,contactInd,posBall_conductive)).^2)));
+
+ C.g_2i = acos(L.B_i/(2*tempR_li));
+ C.g_2a = acos(L.B_a/(2*tempR_la));
+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);
+
+ 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.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? ###
+
+ 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.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;
+ for ii = indices
+ x1=double(C.g_1i(ii));
+ x2=double(C.g_2i);
+
+ C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_i(p)./C.h_i(p,t,ii), x1 , x2 , 0, pi/9);
+ C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) C.dA_a(p)./C.h_a(p,t,ii), C.g_1a(ii)+pi, C.g_2a+pi, 0, pi/9);
+% Z.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_i(p)./Z.h_i(p,t,ii), Z.g_1i(ii) , Z.g_2i , 0, pi/9);
+% Z.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(p,t) Z.dA_a(p)./Z.h_a(p,t,ii), Z.g_1a(ii)+pi, Z.g_2a+pi, 0, pi/9);
+ end
+end
+
+function tobias_steffen(indices)
+ C.h_i = @(alpha,beta,ii) (sqrt(tempR_li^2 - C.DeltaR_i(ii)^2*sin(alpha).^2) - C.DeltaR_i(ii)*cos(alpha)) ./ cos(beta) - tempR_RE;
+ C.h_a = @(alpha,beta,ii) (sqrt(tempR_la^2 - C.DeltaR_a(ii)^2*sin(alpha).^2) - C.DeltaR_a(ii)*cos(alpha)) ./ cos(beta) - tempR_RE;
+ if strcmp(method.unloadedRE,'TobiasSteffen_Kugelfläche')
+ C.dA_i = @(~,beta,~) tempR_RE^2 * cos(beta);
+ C.dA_a = C.dA_i;
+ else
+ C.dA_i = @(alpha,beta,ii) tempR_li./cos(beta) .* (1 - C.DeltaR_i(ii)*cos(alpha)./sqrt(tempR_li^2-C.DeltaR_i(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_i(alpha,beta,ii));
+ C.dA_a = @(alpha,beta,ii) tempR_la./cos(beta) .* (1 - C.DeltaR_a(ii)*cos(alpha)./sqrt(tempR_la^2-C.DeltaR_a(ii)^2*sin(alpha).^2)) .* (tempR_RE + C.h_a(alpha,beta,ii));
+ end
+ for ii = indices
+ C.C_out(1,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_i(alpha,beta,ii)./C.h_i(alpha,beta,ii), C.phi_1i(ii), C.phi_2i, 0, pi/3);
+ C.C_out(2,ii) = 4 * obj.epsilon_0.*S.epsilon_Oel.*integral2(@(alpha,beta) C.dA_a(alpha,beta,ii)./C.h_a(alpha,beta,ii), C.phi_1a(ii), C.phi_2a, 0, pi/3);
+ end
+end
+
+function runSemianalytisch3D(indices)
+
+ if isfield(C,'C_out')
+ C_semi_i = C.C_out(1,:,posBall_conductive);
+ C_semi_a = C.C_out(2,:,posBall_conductive);
+ else
+ C_semi_i = nan( numel(temp_s)/2, 1);
+ C_semi_a = nan( numel(temp_s)/2, 1);
+ end
+
+ 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);
+
+ parfor runSemi = indices
+ C_semi_i(runSemi) = feval(tmp_semi_i, runSemi); %da parfor loop nicht auf nested functions (semianalytisch3D) direkt zugreifen kann
+ C_semi_a(runSemi )= feval(tmp_semi_a, runSemi);
+ end
+
+ 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);
+ end
+ end
+
+ C.C_out(1,:,posBall_conductive)=C_semi_i;
+ C.C_out(2,:,posBall_conductive)=C_semi_a;
+end
+
+function C_Zusatz=semianalytisch3D(s,R_WK,R_R,B_R,B,R_L,epsilon_r,varargin)
+% Semianalytische Berechnung der Kapazität von unbelasteten Wälzkörpern
+% sowie Randbereichen belasteter WK am Innen und Außenring
+% Autor: Steffen Puchtler
+% Herleitung: Masterthesis "Optimierung des Berechnungsverfahrens für die
+% elektrische Kapazität von EHD-Kontakten unter Berücksichtigung des realen
+% elektrischen Feldes", S. 23-26
+%
+% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r)
+% belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b)
+
+% Diskretisierung nicht Default:
+% unbelastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,0,0,nt,np)
+% belastet: C_out = semianalytisch3D(S,R_WK,f,B_R,B,R_L,epsilon_r,a,b,nt,np)
+
+% Inputparameter:
+% S Minimaler Abstand zwischen Wälzkörper und Laufbahn m
+% R_WK Wälzkörperradius m
+% f Schmiegung (R_R / D_WK) -
+% B_R Rillenbreite m
+% B Breite des Lagers m
+% R_L Radius der Laufbahn m !!! R_L < 0 am Innenring
+% epsilon_r Relative Permittivität des Schmierstoffs -
+% a,b Halbachsen der Hertzschen Flächen m
+% nt,np Anzahl der Diskretisierngen in theta bzw. phi -
+
+% profile on
+
+ switch nargin
+ case 7
+ nt = 60; np = 60;
+ belastet = false;
+ case 9
+ nt = 110; np = 110;
+ a = varargin{1};
+ b = varargin{2};
+ belastet = a~=0 && b~=0;
+ case 11
+ a = varargin{1};
+ b = varargin{2};
+ nt = varargin{3};
+ np = varargin{4};
+ belastet = a~=0 && b~=0;
+ otherwise
+ warning('The number of input arguments is wrong! Please check the input of semianalytisch3D.')
+ end
+
+ DeltaR = R_R-R_WK-s;
+ R_T = R_L - R_R;
+ Dz = R_R-sqrt(R_R^2-B_R^2/4);
+ r_yz = @(phi) (R_WK-R_L+s).*cos(phi) + sqrt((R_L-Dz).^2-(R_WK-R_L+s).^2.*sin(phi).^2)*sign(R_L);
+ if belastet
+ Theta_0 = @(phi) a*sqrt(max(1/R_WK^2-phi.^2/b^2,0));
+ else
+ Theta_0 = @(~) 0;
+ end
+ Theta_1 = @(phi) atan(B_R/2./r_yz(phi));
+ Theta_2 = @(phi) atan(B /2./r_yz(phi));
+ if R_L > 0
+ phi_1 = pi/2;
+ else
+ phi_1 = asin(R_L/(R_L-s-R_WK))*0.99;
+ end
+ p = linspace(0,phi_1,np);
+ P = repmat(p',1,nt);
+ T_t = zeros(size(P));
+ for ii = 1:np
+ T_t(ii,:) = linspace(Theta_0(p(ii)),Theta_1(p(ii)),nt);
+ end
+
+ r = zeros(size(T_t));
+ for ii = 1:size(T_t,1)
+ for jj = 1:size(T_t,2)
+ theta = T_t(ii,jj);
+ phi = P(ii,jj);
+ fun = @(r) (R_T+R_R*sqrt(1-(r.*sin(theta)/R_R).^2)).^2 - (r.*sin(phi).*cos(theta)).^2 - (DeltaR+R_T+r.*cos(phi).*cos(theta)).^2;
+% options=optimset('Display','off','FunValCheck','off','OutputFcn',[],'PlotFcns',[],'TolX',eps);
+ r(ii,jj) = fzero(fun,R_WK);
+ end
+ end
+ clear phi theta
+ invalid = r(:)>1;
+ r2small = r(:)/R_WK<1.00125;
+ if any(invalid)
+ warning('%i Entries of r have been removed',sum(invalid))
+ C_Zusatz=nan;
+ return
+ 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);
+ firstInteg = zeros(1,np);
+ for ii = 1:np
+ firstInteg(ii) = trapz(T_t(ii,:),integrand_t(ii,:));
+ end
+ C_Rille = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
+
+
+ T_k = zeros(size(P));
+ for ii = 1:np
+ T_k(ii,:) = linspace(Theta_1(p(ii)),Theta_2(p(ii)),nt);
+ end
+ integrand_k = R_WK^2./((r_yz(P)-s)./cos(T_k)-R_WK).*cos(T_k);
+ firstInteg = zeros(1,np);
+ for ii = 1:np
+ firstInteg(ii) = trapz(T_k(ii,:),integrand_k(ii,:));
+ end
+ C_Rand = 4*epsilon_r * obj.epsilon_0 * trapz(P(:,1)',firstInteg);
+
+ C_Zusatz=C_Rille + C_Rand;
+% if nargout == 1
+% varargout = {C_Rille + C_Rand};
+% elseif nargout == 2
+% varargout = {C_Rille,C_Rand};
+% end
+
+% profile off
+% profile viewer
+
+end
+
+function seperateVecs=splitVec(vec,status)
+ %da die parfor-Schleife nur Vektoren bearbeiten kann, die in Ganzzahl
+ %Schritten auf- oder absteigen, wird in dieser Funktion jeder Vektor in
+ %mehrere Vektoren aufgeteilt, die diese Bedingungen erfüllen.
+
+ findInds = find(vec==status);
+ startInds = unique([findInds(1) findInds(diff([0 findInds])>1)]); % Segment Start Indices
+
+ if length(startInds)==1
+ seperateVecs{1}=findInds;
+ else
+ findNonZero = [find(diff([0 findInds])>1)-1 length(findInds)];
+ findNonZero(findNonZero == 0) = [];
+ endInds = findInds(findNonZero); % Segment End Indices
+
+ for runSection = 1:length(startInds)
+ seperateVecs{runSection} = (startInds(runSection):endInds(runSection));
+ end
+ end
+end
+
+function ind_psi = check_convert_psi_calc2matrx(psi,Z,symCase)
+ % Diese Funktion überprüft, ob der Input-Vektor psi alle Winkel
+ % enthält, die benötigt werden, um die Gesamtkapazität des Wälzlagers
+ % zu berechnen.
+ % Außerdem wird eine Matrix ind_psi erstellt, die die benötigten
+ % Indizes ausgibt, um dieberechneten Einzelkapazitäten in die richtige
+ % Reihenfolge und Ordnung zu bringen.
+
+ switch symCase
+ case 'useSymmetry'
+ symmetry = true;
+ case 'noSymmetry'
+ symmetry = false;
+ otherwise
+ error('Wrong entry for symCase! Possible Input: ''useSymmetry'' or ''noSymmetry''')
+ end
+ ind_psi = nan(obj.L.Z,length(obj.psi));
+
+ psi_calc = psi;
+
+ for run_C = 1 : length(obj.psi)
+
+ psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+ psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+ while max(psiNeedForIndCapa) > (min(obj.psi) + 2*pi - 1e-8)
+ psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) = psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+ end
+
+ if ~isnan(obj.psi_calc)
+
+ psi_check = psiNeedForIndCapa;
+ symPoint = ceil(min(psi)/pi)*pi+pi;
+ if symmetry %Es entsteht eine Symmetrie in der vertikalen Achse, wodurch die Berechnung von zB. der Kapazität nur für eine Hälfte der Umdrehung berechnet werden muss
+ psi_check(psi_check > symPoint) = 2*symPoint - psi_check(psi_check > symPoint);
+ tempPsiComp = 2*symPoint - psi_check(psi_check < symPoint - eps(1e8)); % um Rechenungenauigkeiten durch irrationale Zahl zu vermeiden, wird der Filten erweitert
+
+ psi_check = uniquetol(psi_check, 1e-8);
+ psi_check(ismembertol(psi_check,tempPsiComp,1e-8)) = []; %gespiegelte Punkte löschen
+ end
+ psi_raw = round(obj.psi,8);
+ while (min(psi_check) < min(psi_raw)) || (max(psi_check) >= min(psi_raw) + 2*pi- eps(1e8))
+ psi_check(psi_check<min(psi_raw)) = psi_check(psi_check<min(psi_raw))+ 2*pi;
+ psi_check(psi_check>=min(psi_raw)+2*pi- eps(1e8)) = psi_check(psi_check>=min(psi_raw)+2*pi- eps(1e8))- 2*pi;
+ end
+
+ [~, ind_single_psi] = ismembertol(psi_check, psi_calc,1e-8);
+ [~, ind_comp] = ismembertol(psiNeedForIndCapa, obj.psi_all,1e-8);
+ ind_psi(:,run_C) = obj.ind_psi_all(ind_comp);
+ else
+ [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
+ ind_psi(:,run_C) = ind_single_psi';
+ end
+ %Check, if all nessecary C_el_single are calculated
+ assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!')
+
+ end
+end
+
+function matrx_out = convert_vek2matr(A,B)
+ % Diese Funktion erzeugt aus dem gewünschten Vektor und einer
+ % Index-Matrix eine Matrix, gefüllt mit den Vektoreinträgen in der
+ % Reihenfolge der Matrix-Indizes.
+
+ matrx_out = nan(size(A,1),size(B,2),size(B,1));
+ for run_1A = 1:size(A,1)
+ for run_1B = 1:size(B,1)
+ matrx_out(run_1A,:,run_1B) = A(run_1A,B(run_1B,:));
+ end
+ end
+end
end
\ No newline at end of file
diff --git a/@BearImp/calcClear.m b/@BearImp/calcClear.m
index f0e13550af4abdd02d6b5573c15e2a928563ac08..1a32a7589b9d8f1f3ae3189bae47906e9801f0c6 100644
--- a/@BearImp/calcClear.m
+++ b/@BearImp/calcClear.m
@@ -1,60 +1,60 @@
-function calcClear(obj)
-%2.2 Lagerluft R = R(L,F_r)
-%Berechnet insbesondere das Betriebsspiel des Wälzlagers P_d in Abhängigkeit
-%von u.a. Auslieferungslagerluft, Einbausituation (Passungen),
-%Setzbeträgen, etc.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Julius van der Kuip
-
-%% Parameter prüfen
-assert(obj.up2date.L,'Lager nicht gesetzt')
-L = obj.L; F_r = obj.F_r;
-R = struct;
-
-%% Berechnung
-R.ue_d = 1/3.*L.ei_d + 2/3.*L.es_d - 2/3.*L.EI_d - 1/3.*L.ES_d;
-R.ue_D = 1/3.*L.ei_D + 2/3.*L.es_D - 2/3.*L.EI_D - 1/3.*L.ES_D;
-R.Deltad_ue = R.ue_d.*L.d_i./L.d;
-R.DeltaD_ue = R.ue_D.*L.D ./L.D_G;
-if L.d <= 0.1
- R.gprime = 1e-6;
-else
- R.gprime = 2e-6;
-end
-if L.D <= 0.1
- R.Gprime = 1e-6;
-else
- R.Gprime = 2e-6;
-end
-R.Deltad_W = 0;
-R.DeltaD_W = 0;
-R.Deltad_F = 0.25e-2.*sqrt(L.d./L.B.*F_r.*1e-3).*1e-6;
-R.DeltaD_F = 0.25e-2.*sqrt(L.D./L.B.*F_r.*1e-3).*1e-6;
-R.Deltad_n = 0;
-R.DeltaD_n = 0;
-R.Deltad_T = 0;
-R.DeltaD_T = 0;
-R.DeltaP_W = R.Deltad_ue + R.gprime + R.Deltad_F + R.Deltad_T + R.Deltad_W + R.Deltad_n;
-R.DeltaP_G = R.DeltaD_ue + R.Gprime - R.DeltaD_F + R.DeltaD_T + R.DeltaD_W + R.DeltaD_n;
-R.ue_W = R.ue_d - R.DeltaP_W;
-R.ue_G = R.ue_D - R.DeltaP_G;
-if R.ue_W > 0 % Überprüfe, ob Spiel oder Übermaß vorliegt
- R.sigma_rir = R.ue_W./L.d .* 1./(((1+(L.d./L.d_f).^2)./(1-(L.d./L.d_f).^2)+1./L.ring .nu)./L.ring .E + ...
- ((1+(L.d_i./L.d).^2)./(1-(L.d_i./L.d).^2)-1./L.shaft.nu)./L.shaft.E);
-else
- R.sigma_rir = 0;
-end
-if R.ue_G > 0
- R.sigma_rar = R.ue_G./L.D .* 1./(((1+(L.D_f./L.D).^2)./(1-(L.D_f./L.D).^2)-1./L.ring .nu)./L.ring .E + ...
- ((1+(L.D./L.D_G).^2)./(1-(L.D./L.D_G).^2)+1./L.housing.nu)./L.housing.E);
-else
- R.sigma_rar = 0;
-end
-R.e_ir = 2000 .* L.d./L.ring.E .* (L.d./L.d_f)./(1-(L.d./L.d_f).^2) .* R.sigma_rir .* 1e-3;
-R.e_ar = 2000 .* L.D./L.ring.E .* (L.D_f./L.D)./(1-(L.D_f./L.D).^2) .* R.sigma_rar .* 1e-3;
-
-%% Attribute ändern
-obj.R = R;
-obj.up2date.R = true;
+function calcClear(obj)
+%2.2 Lagerluft R = R(L,F_r)
+%Berechnet insbesondere das Betriebsspiel des Wälzlagers P_d in Abhängigkeit
+%von u.a. Auslieferungslagerluft, Einbausituation (Passungen),
+%Setzbeträgen, etc.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+L = obj.L; F_r = obj.F_r;
+R = struct;
+
+%% Berechnung
+R.ue_d = 1/3.*L.ei_d + 2/3.*L.es_d - 2/3.*L.EI_d - 1/3.*L.ES_d;
+R.ue_D = 1/3.*L.ei_D + 2/3.*L.es_D - 2/3.*L.EI_D - 1/3.*L.ES_D;
+R.Deltad_ue = R.ue_d.*L.d_i./L.d;
+R.DeltaD_ue = R.ue_D.*L.D ./L.D_G;
+if L.d <= 0.1
+ R.gprime = 1e-6;
+else
+ R.gprime = 2e-6;
+end
+if L.D <= 0.1
+ R.Gprime = 1e-6;
+else
+ R.Gprime = 2e-6;
+end
+R.Deltad_W = 0;
+R.DeltaD_W = 0;
+R.Deltad_F = 0.25e-2.*sqrt(L.d./L.B.*F_r.*1e-3).*1e-6;
+R.DeltaD_F = 0.25e-2.*sqrt(L.D./L.B.*F_r.*1e-3).*1e-6;
+R.Deltad_n = 0;
+R.DeltaD_n = 0;
+R.Deltad_T = 0;
+R.DeltaD_T = 0;
+R.DeltaP_W = R.Deltad_ue + R.gprime + R.Deltad_F + R.Deltad_T + R.Deltad_W + R.Deltad_n;
+R.DeltaP_G = R.DeltaD_ue + R.Gprime - R.DeltaD_F + R.DeltaD_T + R.DeltaD_W + R.DeltaD_n;
+R.ue_W = R.ue_d - R.DeltaP_W;
+R.ue_G = R.ue_D - R.DeltaP_G;
+if R.ue_W > 0 % Überprüfe, ob Spiel oder Übermaß vorliegt
+ R.sigma_rir = R.ue_W./L.d .* 1./(((1+(L.d./L.d_f).^2)./(1-(L.d./L.d_f).^2)+1./L.ring .nu)./L.ring .E + ...
+ ((1+(L.d_i./L.d).^2)./(1-(L.d_i./L.d).^2)-1./L.shaft.nu)./L.shaft.E);
+else
+ R.sigma_rir = 0;
+end
+if R.ue_G > 0
+ R.sigma_rar = R.ue_G./L.D .* 1./(((1+(L.D_f./L.D).^2)./(1-(L.D_f./L.D).^2)-1./L.ring .nu)./L.ring .E + ...
+ ((1+(L.D./L.D_G).^2)./(1-(L.D./L.D_G).^2)+1./L.housing.nu)./L.housing.E);
+else
+ R.sigma_rar = 0;
+end
+R.e_ir = 2000 .* L.d./L.ring.E .* (L.d./L.d_f)./(1-(L.d./L.d_f).^2) .* R.sigma_rir .* 1e-3;
+R.e_ar = 2000 .* L.D./L.ring.E .* (L.D_f./L.D)./(1-(L.D_f./L.D).^2) .* R.sigma_rar .* 1e-3;
+
+%% Attribute ändern
+obj.R = R;
+obj.up2date.R = true;
end
\ No newline at end of file
diff --git a/@BearImp/calcFilm.m b/@BearImp/calcFilm.m
index 2a629ef8b38c116566ca3d1931e32887b1054faf..799296c8c09dd24c4ba093f44315afc0c4cb38b3 100644
--- a/@BearImp/calcFilm.m
+++ b/@BearImp/calcFilm.m
@@ -1,157 +1,158 @@
-function varargout = calcFilm(obj, options)
-%2.5 Schmierfilmdicke H = H(S,T,G,B,T_Oil,n)
-%Berechnet die Hertz'schen Flächen sowie die Schmierfilmdicken aller
-%Einzelkontakt (innen und außen für alle belasteten Wälzkörper)
-%
-% options: - 'Q_singleContact' -> ermöglicht die Berechnung der
-% Schmierfilmdicken von externen Lasten
-%
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
-%
-%%%%% Unterschied ein vs kein Ausgabe Argument
-arguments
- obj
- options.Q_inner_singleContact double = obj.B.Q_contInd
- options.Q_outer_singleContact double = obj.B.Q_contInd
- options.BallMaterialInd = NaN % wird nur bei der dynamischen Lastverteilung benötigt
-end
-
-nargoutchk(0,1)
-switch nargout
- case 0
- isPreCalc = false;
- case 1
- isPreCalc = true;
-end
-%% Parameter prüfen
-assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
-assert(obj.up2date.T,'Schmierstoffparameter noch nicht berechnet')
-assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
-if isPreCalc == 0
- assert(obj.up2date.B,'Belastungsverteilung noch nicht berechnet')
-end
-L=obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; T_Oil = obj.T_Oil; omega = obj.omega;
-H = struct;
-
-method = possibleMethods.addDefault(obj.method).H;
-
-%% Berechnung
-
-Q_temp(1,:,:)=options.Q_inner_singleContact';
-Q_temp(2,:,:)=options.Q_outer_singleContact';
-
-
-switch isPreCalc
- case 0
- if L.allBallsSameMaterial
- ballMaterialInd=1;
- else
- ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
- end
- case 1
- ballMaterialInd = options.BallMaterialInd;
-end
-
-H.b = (6*[G.fraktE_i;G.fraktE_a].*Q_temp.*reshape([G.R_si(ballMaterialInd);G.R_sa(ballMaterialInd)],2,1,[]) ./ (pi*[G.k_i;G.k_a].*reshape([G.E_red(ballMaterialInd);G.E_red(ballMaterialInd)],2,1,[]))).^(1/3);
-H.a = [G.k_i;G.k_a].*H.b;
-H.A = pi*H.a.*H.b;
-H.p =Q_temp./H.A;
-H.alpha_p = 1 ./ (S.a_1 + S.a_2*T_Oil + (S.b_1 + S.b_2*T_Oil).*H.p);
-H.u_ir= (G.d_m.^2 - G.D_RE(ballMaterialInd).^2 .* cosd(G.alpha).^2) ./ (4*G.d_m) .* omega;
-H.u_ar = H.u_ir;
-H.u(:,1,:) = [H.u_ir; H.u_ar];
-E_red_tmp(1,1,:) = G.E_red(ballMaterialInd);
-R_x_temp(:,1,:) = G.R_x(:,ballMaterialInd);
-R_y_temp(:,1,:) = G.R_y(:,ballMaterialInd);
-H.U = T.eta_0 .* H.u ./ (E_red_tmp .* R_x_temp);
-H.G = H.alpha_p.*E_red_tmp;
-H.W = Q_temp./(E_red_tmp .* R_x_temp.^2);
-
-switch method
- case 'Hamrock/Dowson'
- k_tmp = [G.k_i;G.k_a];
- H.H_0 = 2.69 * H.G.^0.53 .* H.U.^0.67 .* H.W.^-0.067 .* (1-0.61*exp(-0.73.*k_tmp));
- H.H_min = 3.63 * H.G.^0.49 .* H.U.^0.68 .* H.W.^-0.073 .* (1- exp(-0.68.*k_tmp));
- H.h_0raw = H.H_0 .*R_x_temp; % without correction factors
- H.h_min = H.H_min.*R_x_temp;
- H.L = T.eta_0 .* T.alpha_etaT .* H.u.^2 ./ (4*S.lambda);
- H.C_korr = 3.94 ./ (3.94 + H.L.^0.62);
- H.h_0 = H.h_0raw .*H.C_korr;
- H.h_minth = H.h_min.*H.C_korr;
-
- case 'Moes'
-
- H.lamda = R_x_temp./R_y_temp;
- H.L = H.G .* H.U.^0.25;
- H.N = H.W .* H.lamda.^0.5 .* H.U.^(-0.75);
-
- H.C_RI = 145 * (1+0.796 * H.lamda .^(14/15)) .^(-15/7);
- H.C_RP = 1.29 * (1+0.691 * H.lamda) .^(-2/3);
- H.C_EI = 3.18 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-14/25);
- H.C_EP = 1.48 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-7/20);
-
- H.H_RI = H.C_RI .* H.N.^(-2);
- H.H_RP = H.C_RP .* H.L.^(2/3);
- H.H_EI = H.C_EI .* H.N.^(-2/15);
- H.H_EP = H.C_EP .* H.N.^(-1/12) .* H.L.^(3/4);
- H.s = 3/2.* (1+ exp(-1.2 .* H.H_EI ./ H.H_RI));
-
- H.H_0 = ((H.H_RI.^(3/2) + (H.H_EI.^(-4) + 0.1 .* H.lamda.^4).^(-3/8)).^(2*H.s./3) + (H.H_RP .^(-8) + H.H_EP .^(-8) ).^(-H.s ./8)) .^(1./H.s);
- H.h_0 = H.H_0 .*R_x_temp .* H.U.^(0.5);
-end
-
-switch isPreCalc
- case 1
- varargout{1} = H.h_0;
- case 0
- B.s=zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
-
- B.s(1,B.noContactInd') = B.intersectionBall(B.noContactInd)';
- B.s(2,B.noContactInd') = B.intersectionBall(B.noContactInd)';
-
- switch possibleMethods.addDefault(obj.method).B
- case 'static'
- if length(size(H.h_0)) <= 2
- B.s(1,B. contactInd') = B.intersectionBall(B.contactInd)' + H.h_0(1,B.contactInd)';
- B.s(2,B. contactInd') = B.intersectionBall(B.contactInd)' + H.h_0(2,B.contactInd)'; % s beschreibt den Abstand von WK zu Laufbahnen, unter Berücksichtigung des Schmierfilms in der Lastzone
- else
- B.s(1,B. contactInd') = B.intersectionBall(B.contactInd)' + H.h_0(1,B.contactInd');
- B.s(2,B. contactInd') = B.intersectionBall(B.contactInd)' + H.h_0(2,B.contactInd'); % s beschreibt den Abstand von WK zu Laufbahnen, unter Berücksichtigung des Schmierfilms in der Lastzone
- end
- case 'dynamic'
- B.s(1,B. contactInd') = B.intersectionBall(B.contactInd)';
- B.s(2,B. contactInd') = B.intersectionBall(B.contactInd)';
- end
-
- % Überprüfen, ob die Schmierfilmberechnung innerhalb des vorgegeben
- % Bereichs liegt
- M = H.W.*(2.*H.U).^(-3/4);
- switch method
- case 'Hamrock/Dowson'
- rangeLoadParM = [25,500]; % Quelle: Non-Dimensional Groups, Marian 2020, Fig7
- rangeViscParL = [5,15];
-
- if strcmp(possibleMethods.addDefault(obj.method).B,'dynamic')
- warning('It is not recommended to use Hamrock/Dowson as lubricant film calculation when dynamic load calculation is selected, because the loads become very small outside the load range.')
- end
- case 'Moes'
- rangeLoadParM = [5,1000];
- rangeViscParL = [1,25];
- end
-
- if max(rangeLoadParM) < max(M,[],'all') || min(rangeLoadParM) > min(M,[],'all')
- warning('The force for the lubricant film thickness calculation is outside the permissible range, which means that the calculated lubricant film is no longer meaningful!')
- end
-
- if max(rangeViscParL) < max(H.L,[],'all') || min(rangeViscParL) > min(H.L,[],'all')
- warning('The viscosity parameter for the lubricant film thickness calculation is outside the permissible range, which means that the calculated lubricant film is no longer meaningful!')
- end
-
- %% Attribute ändern
- obj.H = H;
- obj.B = B;
- obj.up2date.H = true;
-end
+function varargout = calcFilm(obj, options)
+%2.5 Schmierfilmdicke H = H(S,T,G,B,T_Oil,n)
+%Berechnet die Hertz'schen Flächen sowie die Schmierfilmdicken aller
+%Einzelkontakt (innen und außen für alle belasteten Wälzkörper)
+%
+% options: - 'Q_singleContact' -> ermöglicht die Berechnung der
+% Schmierfilmdicken von externen Lasten
+%
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Tobias Schirra, Julius van der Kuip
+%
+%%%%% Unterschied ein vs kein Ausgabe Argument
+arguments
+ obj
+ options.Q_inner_singleContact double = obj.B.Q_contInd
+ options.Q_outer_singleContact double = obj.B.Q_contInd
+ options.BallMaterialInd = NaN % wird nur bei der dynamischen Lastverteilung benötigt
+end
+
+nargoutchk(0,1)
+switch nargout
+ case 0
+ isPreCalc = false;
+ case 1
+ isPreCalc = true;
+end
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
+assert(obj.up2date.T,'Schmierstoffparameter noch nicht berechnet')
+assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
+if isPreCalc == 0
+ assert(obj.up2date.B,'Belastungsverteilung noch nicht berechnet')
+end
+L=obj.L; S = obj.S; T = obj.T; G = obj.G; B = obj.B; T_Oil = obj.T_Oil; omega = obj.omega;
+H = struct;
+
+method = possibleMethods.addDefault(obj.method).H;
+
+%% Berechnung
+
+Q_temp(1,:,:)=options.Q_inner_singleContact';
+Q_temp(2,:,:)=options.Q_outer_singleContact';
+
+
+switch isPreCalc
+ case 0
+ if L.allBallsSameMaterial
+ ballMaterialInd=1;
+ else
+ ballMaterialInd=L.RE_order_num(L.IndBall_conductive);
+ end
+ case 1
+ ballMaterialInd = options.BallMaterialInd;
+end
+
+H.b = (6*[G.fraktE_i;G.fraktE_a].*Q_temp.*reshape([G.R_si(ballMaterialInd);G.R_sa(ballMaterialInd)],2,1,[]) ./ (pi*[G.k_i;G.k_a].*reshape([G.E_red(ballMaterialInd);G.E_red(ballMaterialInd)],2,1,[]))).^(1/3);
+H.a = [G.k_i;G.k_a].*H.b;
+H.A = pi*H.a.*H.b;
+H.p =Q_temp./H.A;
+H.alpha_p = 1 ./ (S.a_1 + S.a_2*T_Oil + (S.b_1 + S.b_2*T_Oil).*H.p);
+H.u_ir= (G.d_m.^2 - G.D_RE(ballMaterialInd).^2 .* cosd(G.alpha).^2) ./ (4*G.d_m) .* omega;
+H.u_ar = H.u_ir;
+H.u(:,1,:) = [H.u_ir; H.u_ar];
+E_red_tmp(1,1,:) = G.E_red(ballMaterialInd);
+R_x_temp(:,1,:) = G.R_x(:,ballMaterialInd);
+R_y_temp(:,1,:) = G.R_y(:,ballMaterialInd);
+H.U = T.eta_0 .* H.u ./ (E_red_tmp .* R_x_temp);
+H.G = H.alpha_p.*E_red_tmp;
+H.W = Q_temp./(E_red_tmp .* R_x_temp.^2);
+
+switch method
+ case 'Hamrock/Dowson'
+ k_tmp = [G.k_i;G.k_a];
+ H.H_0 = 2.69 * H.G.^0.53 .* H.U.^0.67 .* H.W.^-0.067 .* (1-0.61*exp(-0.73.*k_tmp));
+ H.H_min = 3.63 * H.G.^0.49 .* H.U.^0.68 .* H.W.^-0.073 .* (1- exp(-0.68.*k_tmp));
+ H.h_0raw = H.H_0 .*R_x_temp; % without correction factors
+ H.h_min = H.H_min.*R_x_temp;
+ H.L = T.eta_0 .* T.alpha_etaT .* H.u.^2 ./ (4*S.lambda);
+ H.C_korr = 3.94 ./ (3.94 + H.L.^0.62);
+ H.h_0 = H.h_0raw .*H.C_korr;
+ H.h_minth = H.h_min.*H.C_korr;
+
+ case 'Moes'
+
+ H.lamda = R_x_temp./R_y_temp;
+ H.L = H.G .* H.U.^0.25;
+ H.N = H.W .* H.lamda.^0.5 .* H.U.^(-0.75);
+
+ H.C_RI = 145 * (1+0.796 * H.lamda .^(14/15)) .^(-15/7);
+ H.C_RP = 1.29 * (1+0.691 * H.lamda) .^(-2/3);
+ H.C_EI = 3.18 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-14/25);
+ H.C_EP = 1.48 * (1+0.006 * log( H.lamda) + 0.63 .* H.lamda .^(4/7)) .^(-7/20);
+
+ H.H_RI = H.C_RI .* H.N.^(-2);
+ H.H_RP = H.C_RP .* H.L.^(2/3);
+ H.H_EI = H.C_EI .* H.N.^(-2/15);
+ H.H_EP = H.C_EP .* H.N.^(-1/12) .* H.L.^(3/4);
+ H.s = 3/2.* (1+ exp(-1.2 .* H.H_EI ./ H.H_RI));
+
+ H.H_0 = ((H.H_RI.^(3/2) + (H.H_EI.^(-4) + 0.1 .* H.lamda.^4).^(-3/8)).^(2*H.s./3) + (H.H_RP .^(-8) + H.H_EP .^(-8) ).^(-H.s ./8)) .^(1./H.s);
+ H.h_0 = H.H_0 .*R_x_temp .* H.U.^(0.5);
+end
+
+switch isPreCalc
+ case 1
+ varargout{1} = H.h_0;
+ case 0
+ switch possibleMethods.addDefault(obj.method).B
+ case 'static'
+
+
+ B.s=zeros(2,B.numOfCagePositions,L.numberOfConductiveBalls);
+
+ B.s(1,B.noContactInd') = B.intersectionBall(1,B.noContactInd)';
+ B.s(2,B.noContactInd') = B.intersectionBall(2,B.noContactInd)';
+
+
+
+ if length(size(H.h_0)) <= 2
+ B.s(1,B. contactInd) = B.intersectionBall(1,B.contactInd) + H.h_0(1,B.contactInd);
+ B.s(2,B. contactInd) = B.intersectionBall(1,B.contactInd) + H.h_0(2,B.contactInd); % s beschreibt den Abstand von WK zu Laufbahnen, unter Berücksichtigung des Schmierfilms in der Lastzone
+ else
+ B.s(1,B. contactInd') = B.intersectionBall(1,B.contactInd') + H.h_0(1,B.contactInd');
+ B.s(2,B. contactInd') = B.intersectionBall(1,B.contactInd') + H.h_0(2,B.contactInd'); % s beschreibt den Abstand von WK zu Laufbahnen, unter Berücksichtigung des Schmierfilms in der Lastzone
+ end
+ end
+
+% Überprüfen, ob die Schmierfilmberechnung innerhalb des vorgegeben
+% Bereichs liegt
+M = H.W.*(2.*H.U).^(-3/4);
+switch method
+ case 'Hamrock/Dowson'
+ rangeLoadParM = [25,500]; % Quelle: Non-Dimensional Groups, Marian 2020, Fig7
+ rangeViscParL = [5,15];
+
+ if strcmp(possibleMethods.addDefault(obj.method).B,'dynamic')
+ warning('It is not recommended to use Hamrock/Dowson as lubricant film calculation when dynamic load calculation is selected, because the loads become very small outside the load range.')
+ end
+ case 'Moes'
+ rangeLoadParM = [5,1000];
+ rangeViscParL = [1,25];
+end
+
+if max(rangeLoadParM) < max(M,[],'all') || min(rangeLoadParM) > min(M,[],'all')
+ warning('The force for the lubricant film thickness calculation is outside the permissible range, which means that the calculated lubricant film is no longer meaningful!')
+end
+
+if max(rangeViscParL) < max(H.L,[],'all') || min(rangeViscParL) > min(H.L,[],'all')
+ warning('The viscosity parameter for the lubricant film thickness calculation is outside the permissible range, which means that the calculated lubricant film is no longer meaningful!')
+end
+end
+%% Attribute ändern
+obj.H = H;
+obj.B = B;
+obj.up2date.H = true;
+
end
\ No newline at end of file
diff --git a/@BearImp/calcGeo.m b/@BearImp/calcGeo.m
index eb6de90608118c69bd5776e6dfb91e222e92db6e..e00841f7a9625b679ab628221e010b0c2673ba72 100644
--- a/@BearImp/calcGeo.m
+++ b/@BearImp/calcGeo.m
@@ -1,102 +1,102 @@
-function calcGeo(obj)
-%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
-%Fläche, Winkelpositionen der Wälzkörper, etc.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Julius van der Kuip
-
-%% 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;
-
-%% Berechnung
-G.alpha = 0;
-
-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.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));
-
-if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnet
- 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));
-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.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;
-G.b_li = b_li*(1+L.ring.alpha_t*G.Delta_Ti_norm);
-G.b_la = b_la*(1+L.ring.alpha_t*G.Delta_Ta_norm);
-
-G.P_d = G.D-G.d-2*G.D_RE-G.b_li-G.b_la;
-
-G.B_i=L.B_i*(1+L.ring.alpha_t*G.Delta_Ti_norm);
-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.R_bi = (G.d + G.b_li) / 2; %Minimaler innerer Rillenradius
-G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius
-
-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.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.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_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;
-
-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;
-
-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)));
-
-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);
-
-G.k_i = fzero(k_i_fkt,5);
-G.k_a = fzero(k_a_fkt,5);
-
-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.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.Deltapsi = 2*pi./L.Z;
-G.psi_ = (0:L.Z-1) .* G.Deltapsi;
-
-%% Attribute ändern
-obj.G = G;
-obj.up2date.G = true;
-
+function calcGeo(obj)
+%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
+%Fläche, Winkelpositionen der Wälzkörper, etc.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+%% 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;
+
+%% Berechnung
+G.alpha = 0;
+
+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.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));
+
+if ~L.allBallsSameMaterial % Materialparameter für alle WK Materialien berechnet
+ 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));
+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.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;
+G.b_li = b_li*(1+L.ring.alpha_t*G.Delta_Ti_norm);
+G.b_la = b_la*(1+L.ring.alpha_t*G.Delta_Ta_norm);
+
+G.P_d = G.D-G.d-2*G.D_RE-G.b_li-G.b_la;
+
+G.B_i=L.B_i*(1+L.ring.alpha_t*G.Delta_Ti_norm);
+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.R_bi = (G.d + G.b_li) / 2; %Minimaler innerer Rillenradius
+G.R_ba = (G.D - G.b_la) / 2; %Maximaler äußererRillenradius
+
+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.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.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_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;
+
+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;
+
+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)));
+
+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);
+
+G.k_i = fzero(k_i_fkt,5);
+G.k_a = fzero(k_a_fkt,5);
+
+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.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.Deltapsi = 2*pi./L.Z;
+G.psi_ = (0:L.Z-1) .* G.Deltapsi;
+
+%% Attribute ändern
+obj.G = G;
+obj.up2date.G = true;
+
end
\ No newline at end of file
diff --git a/@BearImp/calcImp.m b/@BearImp/calcImp.m
index 81babf5bca0d321406b759ffb8185ebd1f64106a..3b13b3f740aa6135249319ad99b03ed1fdbad2df 100644
--- a/@BearImp/calcImp.m
+++ b/@BearImp/calcImp.m
@@ -1,72 +1,72 @@
-function calcImp(obj)
-%2.7 Impedanz Z = Z(L,S,T,G,B,H,C)
-%Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
-%und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wir der
-%Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
-%berechnet.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
-
-%% Parameter prüfen
-assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.C,'Kapazitäten noch nicht berechnet')
-L = obj.L; C = obj.C; psi = obj.psi;
-Z = struct;
-%% Berechnung
-
-% Überprüfen, ob alle C_el_single vorhanden
-if ~isnan(obj.psi_calc)
- psi = obj.psi_calc;
-end
-check_all_nessecary_psi(psi,L.Z,L.allBallsSameMaterial)
-
-%Berechnung der approximierten Impedanz
-Z.omega =@(f) 2 * pi * f; %f = frequency of current
-Z.phi_aprx =@(f) atand( -Z.omega(f) .* Z.C_el_aprx .* Z.R_el_aprx);
-Z.Z_el_aprx =@(f) abs(1./(1./ C.R_el_aprx + 1j * Z.omega(f) .* C.C_el_aprx));
-
-%Berechung der exakten Impedanz
-if strcmp(L.cage_electr,'conductor')
- Z.Z_el_single_iext =@(f) 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:));
- Z.Z_el_single_aext =@(f) 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:));
- Z.Z_el_iext =@(f) 1./ sum( 1 ./ Z.Z_el_single_iext(f),3);
- Z.Z_el_aext =@(f) 1./ sum( 1 ./ Z.Z_el_single_aext(f),3);
- Z.Z_el =@(f) Z.Z_el_iext(f) + Z.Z_el_aext(f);
-elseif strcmp(L.cage_electr,'isolator')
- Z.Z_el_single_iext =@(f) squeeze( 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:)) );
- Z.Z_el_single_aext =@(f) squeeze( 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:)) );
- Z.Z_el_REext =@(f) squeeze( Z.Z_el_single_iext(f) + Z.Z_el_single_aext(f) );
- Z.Z_el =@(f) squeeze( sum(Z.Z_el_REext(f),1) );
-
-else
- warning("The electrical properties of the cage are incorrectly defined! No total impedance can be calculated. ['conductor','isolator']")
-end
-%% Attribute ändern
-obj.Z = Z;
-obj.up2date.Z = true;
-
-%% Helper Functions
-function check_all_nessecary_psi(psi,Z,symCase)
-
- if symCase
- psi_calc = [psi, 2*pi - psi];
- else
- psi_calc = psi;
- end
-
- for run_C = 1 : length(obj.psi)
-
- psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
-
- psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) = psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
-
- psiNeedForIndCapa = round(psiNeedForIndCapa,8);
-
- [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
-
- %Check, if all nessecary C_el_single are calculated
- assert(~any(ind_single_psi==0),'Not all required single capacities were calculated in order to calculate the total capacity for a position!')
- end
-end
+function calcImp(obj)
+%2.7 Impedanz Z = Z(L,S,T,G,B,H,C)
+%Berechnet die Kapazität der Hertz'schen Flächen, sowie der Randbereiche
+%und unbelasteter Wälzkörper je nach gewählter Methode. Zusätzlich wir der
+%Widerstand der Hertz'schen Flächen bestimmt und anschließend die Impedanz
+%berechnet.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Tobias Schirra, Leander, Julius van der Kuip
+
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+assert(obj.up2date.C,'Kapazitäten noch nicht berechnet')
+L = obj.L; C = obj.C; psi = obj.psi;
+Z = struct;
+%% Berechnung
+
+% Überprüfen, ob alle C_el_single vorhanden
+if ~isnan(obj.psi_calc)
+ psi = obj.psi_calc;
+end
+check_all_nessecary_psi(psi,L.Z,L.allBallsSameMaterial)
+
+%Berechnung der approximierten Impedanz
+Z.omega =@(f) 2 * pi * f; %f = frequency of current
+Z.phi_aprx =@(f) atand( -Z.omega(f) .* Z.C_el_aprx .* Z.R_el_aprx);
+Z.Z_el_aprx =@(f) abs(1./(1./ C.R_el_aprx + 1j * Z.omega(f) .* C.C_el_aprx));
+
+%Berechung der exakten Impedanz
+if strcmp(L.cage_electr,'conductor')
+ Z.Z_el_single_iext =@(f) 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:));
+ Z.Z_el_single_aext =@(f) 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:));
+ Z.Z_el_iext =@(f) 1./ sum( 1 ./ Z.Z_el_single_iext(f),3);
+ Z.Z_el_aext =@(f) 1./ sum( 1 ./ Z.Z_el_single_aext(f),3);
+ Z.Z_el =@(f) Z.Z_el_iext(f) + Z.Z_el_aext(f);
+elseif strcmp(L.cage_electr,'isolator')
+ Z.Z_el_single_iext =@(f) squeeze( 1./(1./ C.R_el_single(1,:,:) + 1j * 2*pi* f .* C.C_el_single(1,:,:)) );
+ Z.Z_el_single_aext =@(f) squeeze( 1./(1./ C.R_el_single(2,:,:) + 1j * 2*pi* f .* C.C_el_single(2,:,:)) );
+ Z.Z_el_REext =@(f) squeeze( Z.Z_el_single_iext(f) + Z.Z_el_single_aext(f) );
+ Z.Z_el =@(f) squeeze( sum(Z.Z_el_REext(f),1) );
+
+else
+ warning("The electrical properties of the cage are incorrectly defined! No total impedance can be calculated. ['conductor','isolator']")
+end
+%% Attribute ändern
+obj.Z = Z;
+obj.up2date.Z = true;
+
+%% Helper Functions
+function check_all_nessecary_psi(psi,Z,symCase)
+
+ if symCase
+ psi_calc = [psi, 2*pi - psi];
+ else
+ psi_calc = psi;
+ end
+
+ for run_C = 1 : length(obj.psi)
+
+ psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj.psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+
+ psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) = psiNeedForIndCapa(psiNeedForIndCapa >= min(psi_calc) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+
+ psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+ [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-8);
+
+ %Check, if all nessecary C_el_single are calculated
+ assert(~any(ind_single_psi==0),'Not all required single capacities were calculated in order to calculate the total capacity for a position!')
+ end
+end
end
\ No newline at end of file
diff --git a/@BearImp/calcLoad.m b/@BearImp/calcLoad.m
index e738e35d10d5a7856aed55641fe01b92ed3a5924..093be59c1a5ca5d2f52629f8492972556d7d0c4f 100644
--- a/@BearImp/calcLoad.m
+++ b/@BearImp/calcLoad.m
@@ -1,221 +1,325 @@
-function calcLoad(obj)
-%2.4 Belastungsverteilung B = B(L,G,F_r,psi)
-% Berechnet die Belastungsverteilung auf die einzelnen Wälzkörper.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Julius van der Kuip
-
-%% Parameter prüfen
-assert(obj.up2date.L,'Lager nicht gesetzt')
-assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
-L = obj.L; G = obj.G; F_r = obj.F_r; AddOn = obj.AddOn; psi=obj.psi;
-B = struct;
-
-B.method = possibleMethods.addDefault(obj.method).B;
-%% Berechnung
-
-if ~isnan(obj.psi_calc)
- psi = obj.psi_calc; % alle folgenden Rechnungen werden mit dem reduzierten Psi berechnet, um alle benötigten Werte zu berechnen
-end
-
-if strcmp(B.method,'dynamic')
- n_ir = obj.omega * 60 / (2*pi);
- n_ar = 0;
- B.v_cage = pi/60*G.d_m.*((1-cos(G.alpha)./(G.d_m./G.D_RE)).*n_ir./2+(1+cos(G.alpha)./(G.d_m./G.D_RE)).*n_ar./2); %Umfangsgeschwindigkeit des Käfigs [m/s] Brändlein S.88; Dowson u_i/o/c S.58
-
- if isfield(L,'RE_B')
- B.rhoRE = [L.RE_A.rho,L.RE_B.rho];
- else
- B.rhoRE = L.RE_A.rho;
- end
-
- B.m_RE = 4./3.*pi.*(L.D_RE/2).^3.*B.rhoRE;
- B.Q_centri = B.m_RE.*B.v_cage.^2./(G.d_m./2); % Zentripetalkraft
-
- B.Q_intp1 = logspace(log10(1e-15),log10(1000*F_r),360);%15
-
- B.delta = nan(L.Z,length(B.Q_intp1));
- for IndBall = 1 : L.Z
- h_0_Film = obj.calcFilm('Q_inner_singleContact',B.Q_intp1,'Q_outer_singleContact',B.Q_intp1+B.Q_centri(L.RE_order_num(IndBall)),'BallMaterialInd',L.RE_order_num(IndBall));
- delta_RE_contact =@(K_p_single,Q_psi,IndBall) ((Q_psi+B.Q_intp1)/K_p_single(L.RE_order_num(IndBall))).^(2/3);
- B.delta(IndBall,:) = G.D_RE(L.RE_order_num(IndBall)) - delta_RE_contact(G.K_pa,B.Q_centri(L.RE_order_num(IndBall)),IndBall) - delta_RE_contact(G.K_pi,0,IndBall) + sum(h_0_Film,1);
- end
-end
-
-B.numOfCagePositions = length(psi);
-B.delta_s=nan(B.numOfCagePositions,2);
-B.delta_s=deltaShaft(psi,F_r,B,G,L,B.delta_s,AddOn)';
-
-if L.allBallsSameMaterial
-
- for IndPosBearing=1:B.numOfCagePositions
- B.Q(IndPosBearing) = totalReactionForcePerBallForYandX(B,L,G,B.delta_s(:,IndPosBearing),L.IndBall_conductive,psi(IndPosBearing));
- end
-
-else
-
- B.Q=zeros(L.numberOfConductiveBalls,B.numOfCagePositions);
-
- for posBall_conductive=1:L.numberOfConductiveBalls
-
- B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
- (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
- B.extendedTotalReactionForce=true;
-
- for IndPosBearing=1:B.numOfCagePositions
- B.Q(posBall_conductive,IndPosBearing)=totalReactionForcePerBallForYandX(B,L,G,B.delta_s(:,IndPosBearing),L.IndBall_conductive(posBall_conductive),B.psi_conductive(posBall_conductive,IndPosBearing));
- end
- end
-end
-
-B. contactInd = B.Q~=0; % Indizes der Wälzkörper, die belastet sind (Q ungleich 0)
-B.noContactInd = B.Q==0; % Indizes der Wälzkörper, die nicht belastet sind
-B.Q_contInd=B.Q; %##um leichter B.Q_contInd anzupassen, falls doch nur ~=0 berechnet werden sollen
-
- for posBall_conductive=1:L.numberOfConductiveBalls
- B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
- (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
- for IndPosBall=1:B.numOfCagePositions
- B.intersectionBall(posBall_conductive,IndPosBall) = (deflectionForYandX(L,G,L.IndBall_conductive(posBall_conductive),B.delta_s(:,IndPosBall),B.psi_conductive(posBall_conductive,IndPosBall))-G.D_RE(L.RE_order_num(posBall_conductive)))./2;
- end %Beschreibt die Überschneidung der WK mit den Laufbahnen (<0 -> Kraft wird übertragen; >0 -> keine Kraft wird übertragen)
- end
-
-if strcmp(B.method,'dynamic')
- if max(B.Q_intp1) < max(B.Q) || min(B.Q_intp1) > min(B.Q)
- error('The interpolated lubricant film thicknesses are extrapolated outside the calculated range! This leads to inaccurate results.')
- end
-end
-
-clear B.method
-%% Attribute ändern
-obj.B = B;
-obj.up2date.B = true;
-
-end
-
-
-function delta_s=deltaShaft(psi,F_r,B,G,L,delta_s,AddOn)
-
- for I=1:B.numOfCagePositions
- B.psi_run=G.psi_+psi(I);
- delta_s=solveEquilibrOfForces(I,F_r,B,G,L,delta_s,AddOn);
- end
-end
-
-function delta_s=solveEquilibrOfForces(I,F_r,B,G,L,delta_s,AddOn)
-
- start= [0.001;0.0001]; % Achtung: die Werte dürfen nicht identisch sein, da im sonderfall atan(1)=psi,
- %es zu Fehlern kommt! (in deflectionForYandX wird dann 0 durch 0 geteilt!
-
- if AddOn.Optimization_Toolbox
- options = optimoptions('fsolve', 'StepTolerance', 1e-10, 'FunctionTolerance', 1e-10, 'FunValCheck', 'on', 'OptimalityTolerance',1e-10,'Display','off');%, 'SpecifyObjectiveGradient', true );%,'Algorithm', 'levenberg-marquardt' );%,"PlotFcns",@optimplotfval);
- delta_s(I,:) = fsolve(@(delta_s) EoF(delta_s,B,G,L,F_r,AddOn), start,options);
-
- else
-
- options = optimset('MaxFunEvals',10000,"MaxIter", 10000,"TolFun", 1e-8,"TolX", 1e-8);%,"PlotFcns",@optimplotfval);
- delta_s(I,:) = fminsearch(@(delta_s) EoF(delta_s,B,G,L,F_r,AddOn), start,options);
- end
-
-end
-
-function EoF = EoF(delta_s,B,G,L,F_r,AddOn)
-
- if delta_s(2) == 0
-
- delta_s = delta_s(1);
-
- EoF = abs(sum(yReactionForceOnlyY(B,G,L,delta_s))-F_r);
-
- elseif AddOn.Optimization_Toolbox
- EoF= [F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)); ...
- sum(xReactionForcePerBallForYandX(B,G,L,delta_s))];
- else
- EoF = abs(F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)))+ ...
- abs(sum(xReactionForcePerBallForYandX(B,G,L,delta_s)));
- end
-end
-
-
-
-function Q_y_y=yReactionForceOnlyY(B,G,L,delta_s)
- Q_y_y=zeros(1,L.Z);
-
- for IndBall=1:L.Z
-
- Q_total_y_perball_cont = totalReactionForcePerBallOnlyY(G,IndBall,delta_s,B.psi_run);
-
- psi_k_y = B.psi_run(IndBall)+atan(sin(B.psi_run(IndBall))*delta_s/(G.R_ba-G.R_RE)-cos(B.psi_run(IndBall))*delta_s);
-
- Q_y_y(IndBall)=Q_total_y_perball_cont*cos(psi_k_y);
- end
-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?
- Q_total_y_perball_cont = abs((delta_nury<0)*(delta_nury))^(3/2)*G.K_p;
-end
-
-
-
-function Q_y_yandx_perball=yReactionForcePerBallForYandX(B,G,L,delta_s)
- Q_y_yandx_perball=zeros(1,L.Z);
-
- for IndBall=1:L.Z
- Q_y_yandx_perball(IndBall)=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,B.psi_run)* ...
- cos(deflectionAngleForYandX(B,L,G,IndBall,delta_s));
- end
-end
-
-function Q_total_yandx_perball=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,psi_run)
- if length(psi_run)==1
- psi_run(IndBall)=psi_run;
- end
-
- switch B.method
- case 'static'
- delta_defl = deflectionForYandX(L,G,IndBall,delta_s,psi_run(IndBall))-G.D_RE(L.RE_order_num(IndBall));
-
- if delta_defl < 0
- Q_total_yandx_perball=abs(delta_defl).^(3/2).*G.K_p(L.RE_order_num(IndBall));
- else
- Q_total_yandx_perball = 0;
- end
-
- case 'dynamic'
- delta_exact = deflectionForYandX(L,G,IndBall,delta_s,psi_run(IndBall));
-
- Q_total_yandx_perball = interp1(B.delta(IndBall,:),B.Q_intp1,delta_exact,'spline','extrap'); % 'spline' extrapoliert negative Werte für kleine delta
- end
-end
-
-function delta_yandx=deflectionForYandX(L,G,IndBall,delta_s,psi_run)
-
- atanxy = atan(delta_s(2,:)./delta_s(1,:));
-
- delta_yandx = sin(psi_run-atanxy).*delta_s(2,:)./(sin(atanxy).*(sin(atan(sin(psi_run-atanxy).*delta_s(2,:)./...
- sin(atanxy)./(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(psi_run-atanxy) ...
- ./sin(atanxy).*delta_s(2,:))))))-G.R_bi+G.R_RE(L.RE_order_num(IndBall));
-
-end
-
-function psi_k_yandx=deflectionAngleForYandX(B,L,G,IndBall,delta_s)
-
- atanxy = atan(delta_s(2)/delta_s(1));
- psi = B.psi_run(IndBall);
-
- psi_k_yandx = psi+atan(sin(psi-atanxy)*delta_s(2)/...
- sin(atanxy)/(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(psi-atanxy)...
- *delta_s(2)/sin(atanxy))); %(bei radialer und axialer Verschiebung)
-end
-
-
-
-function Q_x_yandx_perball=xReactionForcePerBallForYandX(B,G,L,delta_s)
- Q_x_yandx_perball=zeros(1,L.Z);
-
- for IndBall=1:L.Z
- Q_x_yandx_perball(IndBall)=totalReactionForcePerBallForYandX(B,L,G,delta_s,IndBall,B.psi_run)*sin(deflectionAngleForYandX(B,L,G,IndBall,delta_s));
- end
+function calcLoad(obj)
+%2.4 Belastungsverteilung B = B(L,G,F_r,psi)
+% Berechnet die Belastungsverteilung auf die einzelnen Wälzkörper.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+%% Parameter prüfen
+assert(obj.up2date.L,'Lager nicht gesetzt')
+assert(obj.up2date.G,'Lagergeometrie in Abhängigkeit des Lagerspiels noch nicht berechnet')
+L = obj.L; G = obj.G; F_r = obj.F_r; F_a = obj.F_a; AddOn = obj.AddOn; psi=obj.psi;
+B = struct;
+
+B.method = possibleMethods.addDefault(obj.method).B;
+%% Berechnung
+
+if ~isnan(obj.psi_calc)
+ psi = obj.psi_calc; % all following calculations are done with the reduced Psi to calculate all needed values
+end
+
+% delta_intp1 (the distance between the raceways) needs to be calculated
+% before delta_s in the dynamic case
+if strcmp(B.method,'dynamic')
+ n_ir = obj.omega * 60 / (2*pi);
+ n_ar = 0;
+ B.v_cage = pi/60*G.d_m.*((1-cos(G.alpha)./(G.d_m./G.D_RE)).*n_ir./2+(1+cos(G.alpha)./(G.d_m./G.D_RE)).*n_ar./2); %Umfangsgeschwindigkeit des Käfigs [m/s] Brändlein S.88; Dowson u_i/o/c S.58
+
+ if isfield(L,'RE_B')
+ B.rhoRE = [L.RE_A.rho,L.RE_B.rho];
+ else
+ B.rhoRE = L.RE_A.rho;
+ end
+
+ B.m_RE = 4./3.*pi.*(L.D_RE/2).^3.*B.rhoRE; % mass in kg
+ B.Q_centri = B.m_RE.*B.v_cage.^2./(G.d_m./2); % Zentripetalkraft
+
+ B.Q_intp1 = logspace(log10(1e-31),log10(1000*sqrt(F_r^2+F_a^2)),360);%15
+
+ B.delta_intp1 = nan(L.Z,length(B.Q_intp1));
+ for IndBall = 1 : L.Z
+ h_0_Film = obj.calcFilm('Q_inner_singleContact',B.Q_intp1,'Q_outer_singleContact',B.Q_intp1+B.Q_centri(L.RE_order_num(IndBall)),'BallMaterialInd',L.RE_order_num(IndBall));
+ delta_RE_contact =@(K_p_single,Q_psi,Q_intp1,IndBall) ((Q_psi+Q_intp1)/K_p_single(L.RE_order_num(IndBall))).^(2/3);
+ B.delta_intp1(IndBall,:) = G.D_RE(L.RE_order_num(IndBall)) - delta_RE_contact(G.K_pa,B.Q_centri(L.RE_order_num(IndBall)),B.Q_intp1,IndBall) - delta_RE_contact(G.K_pi,0,B.Q_intp1,IndBall) + sum(h_0_Film,1);
+ end
+end
+
+B.numOfCagePositions = length(psi);
+
+% delta_s (shaft displacement) can be simplified in case the axial force = 0 (than only
+% y- and x- dircetion are of interest)
+if F_a == 0
+ B.delta_s=nan(B.numOfCagePositions,2);
+else
+ B.delta_s=nan(B.numOfCagePositions,3);
+end
+
+B.delta_s=deltaShaft(psi,F_r,F_a,B,G,L,B.delta_s,AddOn)';
+
+% Q (reaction force) needs to be calculated seperatly because it is not saved
+% while calculating delta_s
+if L.allBallsSameMaterial
+
+ for IndPosBearing=1:B.numOfCagePositions
+ B.Q(IndPosBearing) = totalReactionForcePerBallForXYZ(B,L,G,B.delta_s(:,IndPosBearing),L.IndBall_conductive,psi(IndPosBearing));
+ end
+
+else
+
+ B.Q=zeros(L.numberOfConductiveBalls,B.numOfCagePositions);
+
+ for posBall_conductive=1:L.numberOfConductiveBalls
+
+ B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
+ (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
+ B.extendedTotalReactionForce=true;
+
+ for IndPosBearing=1:B.numOfCagePositions
+ B.Q(posBall_conductive,IndPosBearing)=totalReactionForcePerBallForXYZ(B,L,G,B.delta_s(:,IndPosBearing),L.IndBall_conductive(posBall_conductive),B.psi_conductive(posBall_conductive,IndPosBearing));
+ end
+ end
+end
+
+if strcmp(B.method,'dynamic')
+ assert(any(max(B.Q_intp1) > max(B.Q) & min(B.Q_intp1) < min(B.Q)), 'The interpolated lubricant film thicknesses are extrapolated outside the calculated range! This leads to inaccurate results.')
+end
+
+B. contactInd = B.Q~=0; % Indices of rolling elements under load (Q not equal to 0)
+B.noContactInd = B.Q==0; % Indices of rolling elements that are not loaded
+B.Q_contInd=B.Q; %##um leichter B.Q_contInd anzupassen, falls doch nur ~=0 berechnet werden sollen
+
+for posBall_conductive=1:L.numberOfConductiveBalls
+ B.psi_conductive(posBall_conductive,:)=[psi(1+round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z):B.numOfCagePositions) psi(1: ...
+ (round((L.IndBall_conductive(posBall_conductive)-1)*B.numOfCagePositions/L.Z)))];
+
+ % intersectionBall (intersection of Ball and inner/outer racerway)
+ if strcmp(B.method,'static') % static assumption: Ball in middle of gap
+ for IndPosBearing=1:B.numOfCagePositions
+ deltaXYZ = deflectionForXYZ(L,G,L.IndBall_conductive(posBall_conductive),B.delta_s(:,IndPosBearing),B.psi_conductive(posBall_conductive,IndPosBearing));
+ B.intersectionBall(1,IndPosBearing,posBall_conductive) = (deltaXYZ - G.D_RE(L.RE_order_num(posBall_conductive)))./2;
+ B.intersectionBall(2,IndPosBearing,posBall_conductive) = B.intersectionBall(1,IndPosBearing,posBall_conductive);
+ B.alpha(posBall_conductive,IndPosBearing) = deflectionAngleForXYZ(L,G,B.delta_s(:,IndPosBearing),deltaXYZ,posBall_conductive);
+
+ end % Describes the overlap of RE with raceways (<0 -> force is transmitted; >0 -> no force is transmitted).
+
+ else % dynamic: Q_centri changes position of ball between raceways
+ for IndPosBearing=1:B.numOfCagePositions
+ h_film = obj.calcFilm('Q_inner_singleContact',B.Q(posBall_conductive,IndPosBearing),'Q_outer_singleContact',B.Q(posBall_conductive,IndPosBearing)+B.Q_centri(L.RE_order_num(IndBall)),'BallMaterialInd',L.RE_order_num(IndBall));
+ B.intersectionBall(:,IndPosBearing,posBall_conductive) = h_film ...
+ - [delta_RE_contact(G.K_pa,B.Q_centri(L.RE_order_num(IndBall)),B.Q(posBall_conductive,IndPosBearing),IndBall); delta_RE_contact(G.K_pi,0,B.Q(posBall_conductive,IndPosBearing),IndBall)];
+
+ B.s(:,IndPosBearing,posBall_conductive) = h_film + B.intersectionBall(:,IndPosBearing,posBall_conductive);
+
+ deltaXYZ = deflectionForXYZ(L,G,L.IndBall_conductive(posBall_conductive),B.delta_s(:,IndPosBearing),B.psi_conductive(posBall_conductive,IndPosBearing));
+ B.alpha(posBall_conductive,IndPosBearing) = deflectionAngleForXYZ(L,G,B.delta_s(:,IndPosBearing),deltaXYZ,IndBall);
+
+ end
+ end
+
+end
+
+
+clear B.method
+%% Attribute ändern
+obj.B = B;
+obj.up2date.B = true;
+
+end
+
+
+function delta_s=deltaShaft(psi,F_r,F_a,B,G,L,delta_s,AddOn)
+
+ for I=1:B.numOfCagePositions
+ B.psi_run=G.psi_+psi(I);
+ delta_s=solveEquilibrOfForces(I,F_r,F_a,B,G,L,delta_s,AddOn);
+ end
+end
+
+function delta_s=solveEquilibrOfForces(I,F_r,F_a,B,G,L,delta_s,AddOn)
+
+ if size(delta_s,2) == 2 % 1d/2d shaft displacement
+ start= [0.001;0.0001]; % Achtung: die Werte dürfen nicht identisch sein, da im sonderfall atan(1)=psi,
+ %es zu Fehlern kommt! (in deflectionForYX wird dann 0 durch 0 geteilt!
+ else % 3d shaft displacement
+ start= [0.0001;0.001;0.00001];
+ end
+
+ if AddOn.Optimization_Toolbox
+ options = optimoptions('fsolve', 'StepTolerance', 1e-10, 'FunctionTolerance', 1e-10, 'FunValCheck', 'on', 'OptimalityTolerance',1e-10,'Display','off');%, 'SpecifyObjectiveGradient', true );%,'Algorithm', 'levenberg-marquardt' );%,"PlotFcns",@optimplotfval);
+ delta_s(I,:) = fsolve(@(delta_s) EoF(delta_s,B,G,L,F_r,F_a,AddOn), start,options);
+
+ else
+
+ options = optimset('MaxFunEvals',10000,"MaxIter", 10000,"TolFun", 1e-8,"TolX", 1e-8);%,"PlotFcns",@optimplotfval);
+ delta_s(I,:) = fminsearch(@(delta_s) EoF(delta_s,B,G,L,F_r,F_a,AddOn), start,options);
+ end
+
+end
+
+function EoF = EoF(delta_s,B,G,L,F_r,F_a,AddOn)
+
+ % 1d shaft displacement
+ if delta_s(1) == 0 && (size(delta_s,1) < 3 || delta_s(3) == 0)
+
+ delta_s = delta_s(2);
+
+ EoF = abs(sum(yReactionForceOnlyY(B,G,L,delta_s))-F_r);
+
+ % 2d shaft displacement
+ elseif delta_s(1) ~= 0 && (size(delta_s,1) < 3 || delta_s(3) == 0)
+ if AddOn.Optimization_Toolbox
+ EoF = [ sum(xReactionForcePerBallForYandX(B,G,L,delta_s)); ...
+ F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s))];
+ else
+ EoF = abs( sum(xReactionForcePerBallForYandX(B,G,L,delta_s)))+ ...
+ abs(F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)));
+ end
+
+ % 3d shaft displacement
+ else
+ if AddOn.Optimization_Toolbox
+ EoF = [ sum(xReactionForcePerBallForYandX(B,G,L,delta_s)); ...
+ F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)); ...
+ F_a-sum(zReactionForcePerBallForXYZ (B,G,L,delta_s))];
+ else
+ EoF = abs( sum(xReactionForcePerBallForYandX(B,G,L,delta_s)))+ ...
+ abs(F_r-sum(yReactionForcePerBallForYandX(B,G,L,delta_s)))+ ...
+ abs(F_a-sum(zReactionForcePerBallForXYZ (B,G,L,delta_s)));
+ end
+ end
+end
+
+
+% only shaft displacement in y-direction
+function Q_y_y=yReactionForceOnlyY(B,G,L,delta_s)
+ Q_y_y=zeros(1,L.Z);
+
+ for IndBall=1:L.Z
+
+ Q_total_y_perball_cont = totalReactionForcePerBallOnlyY(G,IndBall,delta_s,B.psi_run);
+
+ psi_k_y = B.psi_run(IndBall)+atan(sin(B.psi_run(IndBall))*delta_s/(G.R_ba-G.R_RE)-cos(B.psi_run(IndBall))*delta_s);
+
+ Q_y_y(IndBall)=Q_total_y_perball_cont*cos(psi_k_y);
+ end
+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?
+ Q_total_y_perball_cont = abs((delta_nury<0)*(delta_nury))^(3/2)*G.K_p;
+end
+
+
+% shaft displacement in y- and x-direction
+function Q_y_yandx_perball=yReactionForcePerBallForYandX(B,G,L,delta_s)
+ Q_y_yandx_perball=zeros(1,L.Z);
+
+ for IndBall=1:L.Z
+ if size(delta_s) == 3 % 3d shaft displacement
+ alpha = deflectionAngleForXYZ(L,G,delta_s,deflectionForYX(L,G,IndBall,delta_s,B.psi_run(IndBall)),IndBall);
+ else % 1d/2d shaft displacement
+ alpha = 0;
+ end
+
+ Q_y_yandx_perball(IndBall)=totalReactionForcePerBallForXYZ(B,L,G,delta_s,IndBall,B.psi_run)* ...
+ cos(deflectionAngleForYandX(B,L,G,IndBall,delta_s))*cos(alpha);
+ end
+end
+
+function Q_x_yandx_perball=xReactionForcePerBallForYandX(B,G,L,delta_s)
+ Q_x_yandx_perball=zeros(1,L.Z);
+
+ for IndBall=1:L.Z
+ if size(delta_s) == 3 % 3d shaft displacement
+ alpha = deflectionAngleForXYZ(L,G,delta_s,deflectionForYX(L,G,IndBall,delta_s,B.psi_run(IndBall)),IndBall);
+ else % 1d/2d shaft displacement
+ alpha = 0;
+ end
+
+ Q_x_yandx_perball(IndBall)=totalReactionForcePerBallForXYZ(B,L,G,delta_s,IndBall,B.psi_run)*...
+ sin(deflectionAngleForYandX(B,L,G,IndBall,delta_s))*cos(alpha);
+ end
+end
+
+function [Q_total_xyz_perball,varargout]=totalReactionForcePerBallForXYZ(B,L,G,delta_s,IndBall,psi_run)
+ if length(psi_run)==1
+ psi_run(IndBall)=psi_run;
+ end
+
+ if size(delta_s,1) == 2 % 2d shaft displacement
+ delta_exact = deflectionForYX(L,G,IndBall,delta_s,psi_run(IndBall)); % distance between inner and outer raceway
+ else % 3d shaft displacement
+ delta_exact = deflectionForXYZ (L,G,IndBall,delta_s,psi_run(IndBall));
+ end
+
+ switch B.method
+ case 'static'
+
+ delta_RWRE_defl = delta_exact - G.D_RE(L.RE_order_num(IndBall)); % distance between raceway(RW) and RE
+
+ if delta_RWRE_defl < 0
+ Q_total_xyz_perball=abs(delta_RWRE_defl).^(3/2).*G.K_p(L.RE_order_num(IndBall));
+ else
+ Q_total_xyz_perball = 0;
+ end
+
+ case 'dynamic'
+ Q_total_xyz_perball = interp1(B.delta_intp1(IndBall,:),B.Q_intp1,delta_exact,'spline','extrap'); % 'spline' extrapoliert fälschlicher Weise negative Werte für kleine delta
+ end
+ if nargout == 2
+ varargout{1} = delta_exact;
+ end
+end
+
+function delta_yandx=deflectionForYX(L,G,IndBall,delta_s,psi_run)
+
+ atanxy = atan(delta_s(1,:)./delta_s(2,:));
+
+ delta_yandx = sin(psi_run-atanxy).*delta_s(1,:)./(sin(atanxy).*(sin(atan(sin(psi_run-atanxy).*delta_s(1,:)./...
+ sin(atanxy)./(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(psi_run-atanxy) ...
+ ./sin(atanxy).*delta_s(1,:))))))-G.R_bi+G.R_RE(L.RE_order_num(IndBall));
+
+end
+
+function psi_k_yandx=deflectionAngleForYandX(B,L,G,IndBall,delta_s)
+
+ atanxy = atan(delta_s(1)/delta_s(2));
+ psi = B.psi_run(IndBall);
+
+ psi_k_yandx = psi+atan(sin(psi-atanxy)*delta_s(1)/...
+ sin(atanxy)/(G.R_ba-G.R_RE(L.RE_order_num(IndBall))-cos(psi-atanxy)...
+ *delta_s(1)/sin(atanxy))); %(bei radialer und axialer Verschiebung)
+end
+
+
+
+% shaft displacement in y-, x- and z-direction
+
+function Q_z_xyz_perball=zReactionForcePerBallForXYZ(B,G,L,delta_s)
+ Q_z_xyz_perball=zeros(1,L.Z);
+
+ for IndBall=1:L.Z
+ Q_z_xyz_perball(IndBall)=totalReactionForcePerBallForXYZ(B,L,G,delta_s,IndBall,B.psi_run)* ...
+ sin(deflectionAngleForXYZ(L,G,delta_s,deflectionForYX(L,G,IndBall,delta_s,B.psi_run(IndBall)),IndBall));
+ end
+end
+
+function delta_XYZ=deflectionForXYZ(L,G,IndBall,delta_s,psi_run)
+
+ if size(delta_s,1) ~= 3
+ delta_s(3) = 0;
+ end
+
+ delta_XYZ = G.R_la(L.RE_order_num(IndBall)) + G.R_li(L.RE_order_num(IndBall))...
+ - sqrt(delta_s(3,:)^2 + (G.R_la(L.RE_order_num(IndBall)) + G.R_li(L.RE_order_num(IndBall))...
+ - deflectionForYX(L,G,IndBall,delta_s,psi_run))^2);
+
+end
+
+function alpha_xyz=deflectionAngleForXYZ(L,G,delta_s,delta,IndBall)
+ if size(delta_s,1) ~= 3
+ delta_s(3) = 0;
+ end
+ alpha_xyz = tan(delta_s(3,:) ./ (G.R_li(L.RE_order_num(IndBall)) + G.R_la(L.RE_order_num(IndBall)) - delta));
end
\ No newline at end of file
diff --git a/@BearImp/calcLub.m b/@BearImp/calcLub.m
index 5a99896a9dea7df14d3eb8813910941f9053b04c..c31043dbbb9d2ae36c48428db80ba68328d4d9bf 100644
--- a/@BearImp/calcLub.m
+++ b/@BearImp/calcLub.m
@@ -1,40 +1,40 @@
-function calcLub(obj)
-%2.1 Schmierstoffparameter (S,T_Oil)
-%Berechnet die temperaturabhängigen Schmierstoffparameter
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
-%% Parameter prüfen
-assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
-S = obj.S; T_Oil = obj.T_Oil;
-T = struct;
-
-T.method = possibleMethods.addDefault(obj.method).T;
-
-%% Berechnung
-T.rho = @(theta) S.rho_15.*(1 - S.alpha_rho.*(theta-15));
-T.eta_040 = T.rho(40).*S.nu_40;
-T.eta_0100 = T.rho(100).*S.nu_100;
-T.alpha_etaT = log(T.eta_040/T.eta_0100)/60;
-switch T.method
- case 'linear'
- T.alpha_eta = (T.eta_0100-T.eta_040)/60;
- 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_0 = T.eta_00 + T.alpha_eta.*T_Oil;
- T.nu_38 = T.nu_0 + T.alpha_nu.*38;
- 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_0 = S.K .* exp(S.B ./ (T_Oil + S.C));
- eta_38 = S.K .* exp(S.B ./ (40 + S.C));
- T.nu_38 = eta_38 ./ T.rho(38);
-end
-
-
-%% Attribute ändern
-obj.T = T;
-obj.up2date.T = true;
-end
-
+function calcLub(obj)
+%2.1 Schmierstoffparameter (S,T_Oil)
+%Berechnet die temperaturabhängigen Schmierstoffparameter
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+%% Parameter prüfen
+assert(obj.up2date.S,'Schmierstoff nicht gesetzt')
+S = obj.S; T_Oil = obj.T_Oil;
+T = struct;
+
+T.method = possibleMethods.addDefault(obj.method).T;
+
+%% Berechnung
+T.rho = @(theta) S.rho_15.*(1 - S.alpha_rho.*(theta-15));
+T.eta_040 = T.rho(40).*S.nu_40;
+T.eta_0100 = T.rho(100).*S.nu_100;
+T.alpha_etaT = log(T.eta_040/T.eta_0100)/60;
+switch T.method
+ case 'linear'
+ T.alpha_eta = (T.eta_0100-T.eta_040)/60;
+ 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_0 = T.eta_00 + T.alpha_eta.*T_Oil;
+ T.nu_38 = T.nu_0 + T.alpha_nu.*38;
+ 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_0 = S.K .* exp(S.B ./ (T_Oil + S.C));
+ eta_38 = S.K .* exp(S.B ./ (40 + S.C));
+ T.nu_38 = eta_38 ./ T.rho(38);
+end
+
+
+%% Attribute ändern
+obj.T = T;
+obj.up2date.T = true;
+end
+
diff --git a/@BearImp/calculate.m b/@BearImp/calculate.m
index 680592b1c986107f19038ca5828424b0a655b8a4..c316479e81a00d394af559ac0bbc1a2fb711c664 100644
--- a/@BearImp/calculate.m
+++ b/@BearImp/calculate.m
@@ -1,126 +1,132 @@
-function calculate (obj, forceMode)
-%Berechnet alle Rechnungsabschnitte, die noch nicht berechnet wurden,
-%sofern alle Eingangsgrößen gesetzt sind.
-%
-%forcemode
-% - 'nonForced' -> Bereits berechnete und noch aktuelle größen werden nicht erneut berechnet (default)
-% - 'forced' -> Alle Größen werden erneut errechnet
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- arguments
- obj
- forceMode {mustBeMember(forceMode,{'nonForced','forced'})} = 'nonForced'
- end
-
- if ~isscalar(obj)
- executeAllObj(obj,@calculate)
- return
- end
-
- switch forceMode
- case 'forced'
- force = true;
- case 'nonForced'
- force = false;
- end
-
- if obj.allUp2date && ~force
- warning('All values already calculated')
- return
- end
- assert(obj.ready2calc,'Not all input parameters set yet')
-
- obj.psi_calc = find_nessecary_psi(obj.psi,obj.L.Z,obj.L.allBallsSameMaterial);
- [obj.psi_all,obj.ind_psi_all] = find_all_ind(obj.psi_calc,obj.L.allBallsSameMaterial);
-
- if ~obj.up2date.T || force
- obj.calcLub
- end
- if ~obj.up2date.R || force
- obj.calcClear
- end
- if ~obj.up2date.G || force
- obj.calcGeo
- end
- if ~obj.up2date.B || force
- obj.calcLoad
- end
- if ~obj.up2date.H || force
- obj.calcFilm
- end
- if ~obj.up2date.C || force
- obj.calcCap
- end
- if ~obj.up2date.Z || force
- obj.calcImp
- end
-
- function psi_calc = find_nessecary_psi(psi,Z, symCase) %(angle, number of balls)
- %this function creates a vector with all positions necessary to
- %calculate the capacity of the bearing
-
- switch symCase
- case 'useSymmetry'
- symmetry = true;
- case 'noSymmetry'
- symmetry = false;
- case obj.L.allBallsSameMaterial % hier wird durch die Eigenschaft, ob alle RE das selbe Material besitzen, die Symmetrie-Eigenschaft des Lagers festgelegt
- symmetry = obj.L.allBallsSameMaterial;
- otherwise
- error('Wrong entry for symCase! Possible Input: ''useSymmetry'' or ''noSymmetry''')
- end
-
-
- psi_calc=[];
-
- for indPsi = 1 : length(psi)
- psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + psi(indPsi); % Bsp. Ergebnis: [2,42,82,122,162,202,242,282,322]/360*2*pi (Offset 2 Grad) im Abstand der WK und dem gewünschten Offsetz von psi
- psi_calc = [psi_calc; psiNeedForIndCapa];
- end
-
- while max(psi_calc) > (min(psi_calc) + 2*pi - 1e-8)
- psi_calc(psi_calc >= min(psi) + 2*pi - 1e-8) = psi_calc(psi_calc >= min(psi) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
- end
-
- psi_calc = round(psi_calc,8);
- psi_calc = uniquetol(psi_calc, 1e-8); %doppelte Werte löschen
-
- symPoint = ceil(min(psi)/pi)*pi+pi;
- if symmetry %Es entsteht eine Symmetrie in der vertikalen Achse, wodurch die Berechnung von zB. der Kapazität nur für eine Hälfte der Umdrehung berechnet werden muss
- psi_calc(psi_calc > symPoint) = 2*symPoint - psi_calc(psi_calc > symPoint);
- tempPsiComp = 2*symPoint - psi_calc(psi_calc < symPoint - eps(1e8)); % um Rechenungenauigkeiten durch irrationale Zahl zu vermeiden, wird der Filten erweitert
-
- psi_calc = uniquetol(psi_calc, 1e-8);
- psi_calc(ismembertol(psi_calc,tempPsiComp,1e-8)) = []; %gespiegelte Punkte löschen
- end
-
-
- end
-
- function [psi_all,ind_psi_all] = find_all_ind(psi,symmetry)
-
- if symmetry
- psi_most = [psi,-psi];
- psi_raw = round(obj.psi,8);
- while (min(psi_most) < min(psi_raw)) || (max(psi_most) >= min(psi_raw) + 2*pi - eps(1e8))
- psi_most(psi_most<min(psi_raw)) = psi_most(psi_most<min(psi_raw))+ 2*pi;
- psi_most(psi_most>=min(psi_raw)+2*pi- eps(1e8)) = psi_most(psi_most>=min(psi_raw)+2*pi- eps(1e8))- 2*pi;
- end
- else
- psi_most = psi;
- end
-
- [psi_most,ind_unique,~] = unique(psi_most,'stable');
- [psi_all, ind_sort] = sort(psi_most);
-
- psi_for_ind = [psi, psi];
- psi_for_ind = psi_for_ind(ind_unique);
- psi_for_ind = psi_for_ind(ind_sort);
-
- for pos = 1: length(psi_all)
- [~,~,ind_psi_all(pos)] = intersect(round(psi_for_ind(pos),8),round(psi,8));
- end
- end
+function calculate (obj, forceMode, options)
+%Berechnet alle Rechnungsabschnitte, die noch nicht berechnet wurden,
+%sofern alle Eingangsgrößen gesetzt sind.
+%
+%forcemode
+% - 'nonForced' -> Bereits berechnete und noch aktuelle größen werden nicht erneut berechnet (default)
+% - 'forced' -> Alle Größen werden erneut errechnet
+%
+%options: - 'calcTo' -> to leave out unnecessary calculations you can
+% choose until which calc... the program should calculate
+%
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+ arguments
+ obj
+ forceMode {mustBeMember(forceMode,{'nonForced','forced'})} = 'nonForced'
+ options.calcTo {mustBeMember(options.calcTo,{'calcLub','calcClear', 'calcGeo', 'calcLoad', 'calcFilm', 'calcCap', 'calcImp'})} = 'calcImp'
+ end
+
+ if ~isscalar(obj)
+ executeAllObj(obj,@calculate)
+ return
+ end
+
+ switch forceMode
+ case 'forced'
+ force = true;
+ case 'nonForced'
+ force = false;
+ end
+
+ if obj.allUp2date && ~force
+ warning('All values already calculated')
+ return
+ end
+ assert(obj.ready2calc,'Not all input parameters set yet')
+
+ obj.psi_calc = find_nessecary_psi(obj.psi,obj.L.Z,obj.L.allBallsSameMaterial);
+ [obj.psi_all,obj.ind_psi_all] = find_all_ind(obj.psi_calc,obj.L.allBallsSameMaterial);
+
+ allCalcList = {'calcLub','calcClear', 'calcGeo', 'calcLoad', 'calcFilm', 'calcCap', 'calcImp'};
+
+ if (~obj.up2date.T || force) && any(strcmp(options.calcTo, allCalcList(1:end)))
+ obj.calcLub
+ end
+ if (~obj.up2date.R || force) && any(strcmp(options.calcTo, allCalcList(2:end)))
+ obj.calcClear
+ end
+ if (~obj.up2date.G || force) && any(strcmp(options.calcTo, allCalcList(3:end)))
+ obj.calcGeo
+ end
+ if (~obj.up2date.B || force) && any(strcmp(options.calcTo, allCalcList(4:end)))
+ obj.calcLoad
+ end
+ if (~obj.up2date.H || force) && any(strcmp(options.calcTo, allCalcList(5:end)))
+ obj.calcFilm
+ end
+ if (~obj.up2date.C || force) && any(strcmp(options.calcTo, allCalcList(6:end)))
+ obj.calcCap
+ end
+ if (~obj.up2date.Z || force) && any(strcmp(options.calcTo, allCalcList(end)))
+ obj.calcImp
+ end
+
+ function psi_calc = find_nessecary_psi(psi,Z, symCase) %(angle, number of balls)
+ %this function creates a vector with all positions necessary to
+ %calculate the capacity of the bearing
+
+ switch symCase
+ case 'useSymmetry'
+ symmetry = true;
+ case 'noSymmetry'
+ symmetry = false;
+ case obj.L.allBallsSameMaterial % hier wird durch die Eigenschaft, ob alle RE das selbe Material besitzen, die Symmetrie-Eigenschaft des Lagers festgelegt
+ symmetry = obj.L.allBallsSameMaterial;
+ otherwise
+ error('Wrong entry for symCase! Possible Input: ''useSymmetry'' or ''noSymmetry''')
+ end
+
+
+ psi_calc=[];
+
+ for indPsi = 1 : length(psi)
+ psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + psi(indPsi); % Bsp. Ergebnis: [2,42,82,122,162,202,242,282,322]/360*2*pi (Offset 2 Grad) im Abstand der WK und dem gewünschten Offsetz von psi
+ psi_calc = [psi_calc; psiNeedForIndCapa];
+ end
+
+ while max(psi_calc) > (min(psi_calc) + 2*pi - 1e-8)
+ psi_calc(psi_calc >= min(psi) + 2*pi - 1e-8) = psi_calc(psi_calc >= min(psi) + 2*pi - 1e-8) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch
+ end
+
+ psi_calc = round(psi_calc,8);
+ psi_calc = uniquetol(psi_calc, 1e-8); %doppelte Werte löschen
+
+ symPoint = ceil(min(psi)/pi)*pi+pi;
+ if symmetry %Es entsteht eine Symmetrie in der vertikalen Achse, wodurch die Berechnung von zB. der Kapazität nur für eine Hälfte der Umdrehung berechnet werden muss
+ psi_calc(psi_calc > symPoint) = 2*symPoint - psi_calc(psi_calc > symPoint);
+ tempPsiComp = 2*symPoint - psi_calc(psi_calc < symPoint - eps(1e8)); % um Rechenungenauigkeiten durch irrationale Zahl zu vermeiden, wird der Filten erweitert
+
+ psi_calc = uniquetol(psi_calc, 1e-8);
+ psi_calc(ismembertol(psi_calc,tempPsiComp,1e-8)) = []; %gespiegelte Punkte löschen
+ end
+
+
+ end
+
+ function [psi_all,ind_psi_all] = find_all_ind(psi,symmetry)
+
+ if symmetry
+ psi_most = [psi,-psi];
+ psi_raw = round(obj.psi,8);
+ while (min(psi_most) < min(psi_raw)) || (max(psi_most) >= min(psi_raw) + 2*pi - eps(1e8))
+ psi_most(psi_most<min(psi_raw)) = psi_most(psi_most<min(psi_raw))+ 2*pi;
+ psi_most(psi_most>=min(psi_raw)+2*pi- eps(1e8)) = psi_most(psi_most>=min(psi_raw)+2*pi- eps(1e8))- 2*pi;
+ end
+ else
+ psi_most = psi;
+ end
+
+ [psi_most,ind_unique,~] = unique(psi_most,'stable');
+ [psi_all, ind_sort] = sort(psi_most);
+
+ psi_for_ind = [psi, psi];
+ psi_for_ind = psi_for_ind(ind_unique);
+ psi_for_ind = psi_for_ind(ind_sort);
+
+ for pos = 1: length(psi_all)
+ [~,~,ind_psi_all(pos)] = intersect(round(psi_for_ind(pos),8),round(psi,8));
+ end
+ end
end
\ No newline at end of file
diff --git a/@BearImp/copyToArray.m b/@BearImp/copyToArray.m
index 0c3938f578738de581a3a9dc02a9b27b9894bd6b..d834224a7b4f5cdfcb98dfecdeb23f85f11842d4 100644
--- a/@BearImp/copyToArray.m
+++ b/@BearImp/copyToArray.m
@@ -1,28 +1,28 @@
-function objArray = copyToArray(obj,varargin)
-%Kopiert ein Objekt in eine Matrix gegebender Größe:
-%objArray = obj.copyToArray(n) > n x n Matrix
-%objArray = obj.copyToArray(m,n) > m x n Matrix
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- switch nargin
- case 1
- objArray = copy(obj);
- case 2
- objArray = copy2D(obj,varargin{1},varargin{1});
- case 3
- objArray = copy2D(obj,varargin{1},varargin{2});
- otherwise
- error('Funktion bisher nur für 2D-Matrizen implementiert')
- end
-end
-
-function objArray = copy2D(obj, m, n)
- objArray(m,n) = BearImp; % preallocate
- for ii = 1:m
- for jj = 1:n
- objArray(ii,jj) = copy(obj);
- end
- end
+function objArray = copyToArray(obj,varargin)
+%Kopiert ein Objekt in eine Matrix gegebender Größe:
+%objArray = obj.copyToArray(n) > n x n Matrix
+%objArray = obj.copyToArray(m,n) > m x n Matrix
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+ switch nargin
+ case 1
+ objArray = copy(obj);
+ case 2
+ objArray = copy2D(obj,varargin{1},varargin{1});
+ case 3
+ objArray = copy2D(obj,varargin{1},varargin{2});
+ otherwise
+ error('Funktion bisher nur für 2D-Matrizen implementiert')
+ end
+end
+
+function objArray = copy2D(obj, m, n)
+ objArray(m,n) = BearImp; % preallocate
+ for ii = 1:m
+ for jj = 1:n
+ objArray(ii,jj) = copy(obj);
+ end
+ end
end
\ No newline at end of file
diff --git a/@BearImp/get.m b/@BearImp/get.m
index 8534fd1d3fc6d1bb2a6e0dd205a30ebaf54e056e..cc0a8bb41b0408ff5f0fd861c946103d7d4496ba 100644
--- a/@BearImp/get.m
+++ b/@BearImp/get.m
@@ -1,70 +1,70 @@
-function value = get(obj, name)
-%Gibt ein Attribut eines BearImp-Objekts aus. Funktioniert
-%auch mit BearImp-Arrays bzw. -Matrizen
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- arguments
- obj BearImp
- name char
- end
-
- if isscalar(obj)
- value = getSingleValue(obj,name);
- else
- value = reproduceValues(obj,name); % preallocate
- for ii = 2:numel(obj)
- if iscell(value)
- value{ii} = getSingleValue(obj(ii),name);
- else
- value(ii) = getSingleValue(obj(ii),name);
- end
- end
- end
-end
-
-function value = getSingleValue(obj, name)
- if issplitable(name)
- value = getSingleValueFromStruct(obj, name);
- else
- value = obj.(name);
- end
-end
-
-function outBool = issplitable(name)
- outBool = numel(strsplit(name,'.')) == 2;
-end
-
-function value = getSingleValueFromStruct(obj, name)
- [property, field] = splitName(name);
- value = obj.(property).(field);
-end
-
-function [property, field] = splitName(name)
- splitCell = strsplit(name,'.');
- switch numel(splitCell)
- case 1
- property = splitCell{1};
- field = '';
- case 2
- property = splitCell{1};
- field = splitCell{2};
- otherwise
- error('Fehler beim Auslesen eines Structes. Es kann maximal auf die erste Unterebene zugegriffen werden, z.B. Z.C_el')
- end
-end
-
-function values = reproduceValues(obj,name)
- sampleValue = getSingleValue(obj(1,1),name);
- if isscalar(sampleValue)
- values = repmat(sampleValue,size(obj));
- else
- values = cell(size(obj));
- values{1} = sampleValue;
- end
-end
-
-function prop = getPropertyFromName(name)
- [prop,~] = splitName(name);
+function value = get(obj, name)
+%Gibt ein Attribut eines BearImp-Objekts aus. Funktioniert
+%auch mit BearImp-Arrays bzw. -Matrizen
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+ arguments
+ obj BearImp
+ name char
+ end
+
+ if isscalar(obj)
+ value = getSingleValue(obj,name);
+ else
+ value = reproduceValues(obj,name); % preallocate
+ for ii = 2:numel(obj)
+ if iscell(value)
+ value{ii} = getSingleValue(obj(ii),name);
+ else
+ value(ii) = getSingleValue(obj(ii),name);
+ end
+ end
+ end
+end
+
+function value = getSingleValue(obj, name)
+ if issplitable(name)
+ value = getSingleValueFromStruct(obj, name);
+ else
+ value = obj.(name);
+ end
+end
+
+function outBool = issplitable(name)
+ outBool = numel(strsplit(name,'.')) == 2;
+end
+
+function value = getSingleValueFromStruct(obj, name)
+ [property, field] = splitName(name);
+ value = obj.(property).(field);
+end
+
+function [property, field] = splitName(name)
+ splitCell = strsplit(name,'.');
+ switch numel(splitCell)
+ case 1
+ property = splitCell{1};
+ field = '';
+ case 2
+ property = splitCell{1};
+ field = splitCell{2};
+ otherwise
+ error('Fehler beim Auslesen eines Structes. Es kann maximal auf die erste Unterebene zugegriffen werden, z.B. Z.C_el')
+ end
+end
+
+function values = reproduceValues(obj,name)
+ sampleValue = getSingleValue(obj(1,1),name);
+ if isscalar(sampleValue)
+ values = repmat(sampleValue,size(obj));
+ else
+ values = cell(size(obj));
+ values{1} = sampleValue;
+ end
+end
+
+function prop = getPropertyFromName(name)
+ [prop,~] = splitName(name);
end
\ No newline at end of file
diff --git a/@BearImp/plot.m b/@BearImp/plot.m
index 38c5ac45eaccad0319960330ff5a300112f31e72..f6fa0708b500fa3915f9ecc669d947dd36f817fc 100644
--- a/@BearImp/plot.m
+++ b/@BearImp/plot.m
@@ -1,247 +1,261 @@
-function plot(obj,name,options)
-% plot -> plottet die wichtigesten Verläufe in einem Fenster
-% plot(name) -> plottet den Verlauf der Variable mit 'name'
-% (mehrere Variablen sind auch möglich mit: name == {'var1','var2',...}
-% Die Namen sind in PlotLabel.csv einsehbar und ergänzbar.)
-%
-% options: - 'plotStyle' -> 'linear'/'polar'
-% - 'figureName' -> name of figure window
-% - 'frequency' -> frequency in HZ (to calculate the impedance a frequency is necessary!)
-% - 'create' -> 'new'/'current' (choose between create new figure or take the already existing one)
-%
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Julius van der Kuip
-
- arguments
- obj
- name = 'Default'
- options.plotStyle char {mustBeMember(options.plotStyle,{'linear','polar'})} = 'linear'
- options.frequency = 2e4
- options.figureName char = 'Ergebnisse BearImp';
- options.create char {mustBeMember(options.create,{'new','current'})} = 'new'
- end
-
- if numel(obj) > 1
- for ii = 1:numel(obj) % TODO: not working with matrices
- obj(ii).plot(name,'plotStyle',options.plotStyle,'create',options.create);
- hold on
- end
- hold off
- return
- end
-
- LabelTable = readtable(obj.PlotLabelTableFile); % Einlesen der PlotLabelTable
-
- if strcmp(name,'Default') % Definiert die Variablen, die bei Default geplottet werden sollen
- var2plot = {'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'};
- elseif ischar(name) % damit spätere Rechnungen veranfacht werden können, wird der char in ein cell verpackt
- var2plot = {name};
- else
- var2plot = name;
- end
-
- % hier wird analysiert, welche Variable geplottet werden soll
- var_pos = nan(1,length(var2plot));
- for run_var = 1: length(var2plot)
- isInTable = strcmp(LabelTable.var_name,var2plot{run_var});
-
- assert(any(find(isInTable,true)),'The entered variable name does not exist. Please check variable name!')
-
- var_pos(run_var) = find(isInTable,true); % var_pos gibt die Indizes an, bei dem die zu plottende Variable abgespeichert ist
- end
-
-
- iter_var = LabelTable.var_name(var_pos);
- len_iter = length(iter_var);
-
- switch options.create
- case 'new'
- figure('name',options.figureName)
- case 'current'
- gcf;
- end
-
- for ii = 1:len_iter
-
- if len_iter > 1 % Erzeugt Subplots, falls mehrere Plots gleichzeitig geplottet werden sollen
- % Optimierung der Fensteranordnung
- l1=(ceil(sqrt(len_iter)))*(ceil(sqrt(len_iter))-1); %Grid 0/-1
- l2=(ceil(sqrt(len_iter))+1)*(ceil(sqrt(len_iter))-1); %Grid nicht Quadrat,sondern -1/+1
-
- if l1 >= len_iter
- subpltx =ceil(sqrt(len_iter));
- subplty = ceil(sqrt(len_iter))-1;
- elseif l2 < len_iter
- l3=ceil(sqrt(len_iter)); %Grid quadratisch 0/0
- subpltx = l3;
- subplty = l3;
- else
- subpltx = ceil(sqrt(len_iter))+1;
- subplty = ceil(sqrt(len_iter))-1;
- end
- % Erzeugung des Subplots
- subplot(subpltx,subplty,ii)
- end
-
- run_var = iter_var(ii);
- row_var = LabelTable(strcmp(LabelTable.var_name,run_var),:); % da nur die ausgewählt Variable bearbeitet wird, wird auch aus der Tabelle nur die entsprechende Reihe ausgewählt
-
- assert(isnumeric(row_var.unit_factor), ['Der ''unit_factor'' der Variable ',row_var.var_name{1}, ' ist keine interpretierbare Zahl'])
- assert(~isnan(row_var.unit_factor), ['Der ''unit_factor'' der Variable ',row_var.var_name{1}, ' ist keine interpretierbare Zahl'])
-
- if strcmp(run_var{1},'Z_el') % da die Impedanz eine Frequenz als Input benötigt, um berechnet zu werden
- assert(~isnan(options.frequency),'To be able to calculate the impedance of the bearing, a frequency must be specified in Hz! \n e.g. '' frequency'', 1e4')
- var_raw = eval(row_var.var{1});
- plot_var_raw = abs(var_raw(options.frequency)) *row_var.unit_factor;
- else
- plot_var_raw = eval(row_var.var{1})*row_var.unit_factor; % hier wird der unit_factor mit berücksichtigt
- end
-
- if size(plot_var_raw,3)>1 % Die Graphen werden immer auf den ersten Wälzkörper bezogen geplottet
- plot_var_raw = squeeze(plot_var_raw(:,:,1));
- end
-
-
- if isnan(obj.psi_calc)
-
- [plot_psi,ind] = check_psi(obj.psi,options.plotStyle);
- plot_var = plot_var_raw(:,ind);
-
- else
-
- % da die Funktion calculate keine Werte doppelt berechnet, wird hier der Vetor wieder vervollständigt
- if strcmp(row_var.var_name,'delta_sx')
- [plot_psi,ind] = check_psi(obj.psi_all,options.plotStyle);
- plot_var = plot_var_raw(:,ind);
- % da delta_sx die einzige Variable ist, die punktsymmetrisch gekürzt wird
-
- plot_var((obj.psi_calc(ind)>pi) & (plot_psi>0)) = -plot_var((obj.psi_calc(ind)>pi) & (plot_psi>0));
- plot_var((obj.psi_calc(ind)<pi) & (plot_psi<0)) = -plot_var((obj.psi_calc(ind)<pi) & (plot_psi<0));
- elseif any(strcmp(row_var.var_name,{'C_el_single','R_el_single','C_el','R_el','Z_el'})) && obj.L.allBallsSameMaterial
- [plot_psi,ind] = check_psi(obj.psi,options.plotStyle);
- plot_var = plot_var_raw(:,ind);
- else
- [plot_psi,ind] = check_psi(obj.psi_all,options.plotStyle);
- plot_var = plot_var_raw(:,ind);
- end
-
- end
-
- switch options.plotStyle
- case 'linear'
- plotLinear(plot_var,plot_psi, row_var.title_var, {row_var.full_name{1};strcat(row_var.var_label{1}," in ",row_var.unit{1})})
- case 'polar'
- plotPolar(plot_var,plot_psi, row_var.title_var, row_var.unit{1})
- end
- end
-
-function plotLinear(vek,psi, pltTitle, pltYlabel)
-% plots a linear graph with all necessary adjustments
- vek(isnan(vek)) = 0;
- plot(psi,vek,'linewidth',1) %length(a.B.delta_s(1,:)
- title(pltTitle)
- xlabel("Cage angle in °")
- ylabel(pltYlabel)
- set(gca,'XTick',(-pi:2*pi/12:pi));
- set(gca,'XTickLabel',(-180:30:181));
- xlim([-pi pi])
- hold off
- grid on
- if size(vek,1)==2
- legend('inner contact','outer contact') %,'Position',[0.76 0.83 0.05 0.07]
- end
-end
-
-function plotPolar(vek,ind, pltTitle, unit)
-% plots a polar graph with all necessary adjustments
- vek_plot = [vek,vek(:,1)];
- vek_plot(isnan(vek_plot)) = 0;
- vek_max = max(vek_plot(1,:));
-
- if vek_max == inf
- temp_vek = vek_plot(~isinf(vek_plot) );
- vek_max = max(temp_vek);
- end
-
- vek_min = min(vek_plot(1,:));
-
- digitsRound = floor(log10(abs(vek_max)))+2;
-
- if abs(digitsRound) == inf
- vek_ax_max = vek_max;
- vek_ax_min = vek_min;
- else
- vek_ax_max = round(vek_max,digitsRound,'significant');
- vek_ax_min = round(vek_min,digitsRound,'significant');
- end
-
-
- if vek_ax_min>vek_min
- vek_diff=round(vek_ax_max-vek_ax_min);
- vek_ax_min = vek_ax_min-vek_diff/5;
- end
-
- if (vek_ax_max-vek_ax_min==0) && (vek_max~=0)
- vek_diff = vek_max*1.1;
- elseif (vek_ax_max-vek_ax_min==0) && (vek_max==0)
- vek_diff = vek_max+0.1;
- else
- vek_diff = (vek_ax_max-vek_ax_min);
- end
-
-
- psi = [ind,ind(1)];
-
- string1 = num2str(linspace(0,360-360/obj.L.Z,obj.L.Z)');
- string2 = nan(obj.L.Z,2);
- for run_RE = 1:obj.L.Z
- string2(run_RE,:) = ' °';
- end
- thetaTickLabel = [string1,string2];
-
- polarplot(psi,vek_plot);
- ax = gca;
- ax.ThetaDir = 'clockwise';
- ax.RAxis.Color = 'black';
- ax.ThetaZeroLocation = 'bottom';
- ax.RLim = [vek_ax_min-vek_diff*5/6 vek_max+vek_diff/4];
- ax.RTick = [vek_ax_min, vek_ax_min+vek_diff/4, vek_ax_min+vek_diff/2, vek_ax_min+vek_diff*3/4, vek_ax_min+vek_diff];
- ax.RTickLabel = {strcat(num2str(vek_ax_min)," ",unit), strcat(num2str(vek_ax_min+vek_diff/4)," ",unit), strcat(num2str(vek_ax_min+vek_diff/2)," ",unit), strcat(num2str(vek_ax_min+vek_diff*3/4)," ",unit), strcat(num2str(vek_ax_max)," ",unit)};
- ax.ThetaTick = linspace(0,360-360/obj.L.Z,obj.L.Z);
- ax.ThetaTickLabel = {thetaTickLabel};
- ax.RAxisLocation = 180;
-
- title(pltTitle)
- hold on
- polarplot(psi,ones(1,length(psi))*vek_ax_min,'black');
- if size(vek,1)==2
- legend('inner contact','outer contact','') %,'Position',[0.83 0.07 0.05 0.07]
- end
- hold off
-end
-
-function [psi_checked,ind] = check_psi(psi,plotStyle)
-
- psi(psi<0) = 2*pi + psi(psi<0);
- psi_raw = psi;
- psi = sort(psi);
-
- ind = nan(1,length(psi));
- for pos = 1: length(psi)
- [~,~,ind(pos)] = intersect(psi(pos),psi_raw);
- end
-
- ind = [ind(psi>=0),ind(psi<0)];
- if ~any(isnan(obj.psi_calc)) && ~(any(strcmp(row_var.var_name,{'C_el_single','R_el_single','C_el','R_el','Z_el'})) && obj.L.allBallsSameMaterial)
- ind = obj.ind_psi_all(ind);
- end
-
- switch plotStyle
- case 'linear'
- psi_checked = [psi(psi>=pi)-2*pi,psi(psi<pi)];
- ind = [ind(psi>=pi),ind(psi<pi)];
- case 'polar'
- psi_checked = psi;
- end
-end
+function plot(obj,name,options)
+% plot -> plottet die wichtigesten Verläufe in einem Fenster
+% plot(name) -> plottet den Verlauf der Variable mit 'name'
+% (mehrere Variablen sind auch möglich mit: name == {'var1','var2',...}
+% Die Namen sind in PlotLabel.csv einsehbar und ergänzbar.)
+%
+% options: - 'plotStyle' -> 'linear'/'polar'
+% - 'figureName' -> name of figure window
+% - 'frequency' -> frequency in HZ (to calculate the impedance a frequency is necessary!)
+% - 'create' -> 'new'/'current'/'new_holdon' (choose between create new figure / take the already existing one / when plotting BearImp Array: first plot creates new figure, all following take current)
+%
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Julius van der Kuip
+
+ arguments
+ obj
+ name = 'Default'
+ options.plotStyle char {mustBeMember(options.plotStyle,{'linear','polar'})} = 'linear'
+ options.frequency = 2e4
+ options.figureName char = 'Ergebnisse BearImp';
+ options.create char {mustBeMember(options.create,{'new','current','new_holdon'})} = 'new'
+ % if you add new option, please remember to insert it to following if-statement
+ end
+
+ if numel(obj) > 1
+
+ if strcmp(options.create, 'new_holdon')
+ obj(1).plot(name,'plotStyle',options.plotStyle,'frequency',options.frequency,...
+ 'figureName',options.figureName,'create','new');
+ hold on
+ startFor = 2;
+ else
+ startFor = 1;
+ end
+
+ for ii = startFor:numel(obj) % TODO: not working with matrices
+ obj(ii).plot(name,'plotStyle',options.plotStyle,'frequency',options.frequency,...
+ 'figureName',options.figureName,'create','current');
+ hold on
+ end
+ hold off
+ return
+ end
+
+ LabelTable = readtable(obj.PlotLabelTableFile); % Einlesen der PlotLabelTable
+
+ if strcmp(name,'Default') % Definiert die Variablen, die bei Default geplottet werden sollen
+ var2plot = {'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'};
+ elseif ischar(name) % damit spätere Rechnungen veranfacht werden können, wird der char in ein cell verpackt
+ var2plot = {name};
+ else
+ var2plot = name;
+ end
+
+ % hier wird analysiert, welche Variable geplottet werden soll
+ var_pos = nan(1,length(var2plot));
+ for run_var = 1: length(var2plot)
+ isInTable = strcmp(LabelTable.var_name,var2plot{run_var});
+
+ assert(any(find(isInTable,true)),'The entered variable name does not exist. Please check variable name!')
+
+ var_pos(run_var) = find(isInTable,true); % var_pos gibt die Indizes an, bei dem die zu plottende Variable abgespeichert ist
+ end
+
+
+ iter_var = LabelTable.var_name(var_pos);
+ len_iter = length(iter_var);
+
+ switch options.create
+ case 'new'
+ figure('name',options.figureName)
+ case 'current'
+ gcf;
+ case 'new_holdon'
+ error('You selected "new_holdon" while plotting only one graph. Create a BearImp Array to use this feature or consider "create" / "new" as option.')
+ end
+
+ for ii = 1:len_iter
+
+ if len_iter > 1 % Erzeugt Subplots, falls mehrere Plots gleichzeitig geplottet werden sollen
+ % Optimierung der Fensteranordnung
+ l1=(ceil(sqrt(len_iter)))*(ceil(sqrt(len_iter))-1); %Grid 0/-1
+ l2=(ceil(sqrt(len_iter))+1)*(ceil(sqrt(len_iter))-1); %Grid nicht Quadrat,sondern -1/+1
+
+ if l1 >= len_iter
+ subpltx =ceil(sqrt(len_iter));
+ subplty = ceil(sqrt(len_iter))-1;
+ elseif l2 < len_iter
+ l3=ceil(sqrt(len_iter)); %Grid quadratisch 0/0
+ subpltx = l3;
+ subplty = l3;
+ else
+ subpltx = ceil(sqrt(len_iter))+1;
+ subplty = ceil(sqrt(len_iter))-1;
+ end
+ % Erzeugung des Subplots
+ subplot(subpltx,subplty,ii)
+ end
+
+ run_var = iter_var(ii);
+ row_var = LabelTable(strcmp(LabelTable.var_name,run_var),:); % da nur die ausgewählt Variable bearbeitet wird, wird auch aus der Tabelle nur die entsprechende Reihe ausgewählt
+
+ assert(isnumeric(row_var.unit_factor), ['Der ''unit_factor'' der Variable ',row_var.var_name{1}, ' ist keine interpretierbare Zahl'])
+ assert( ~isnan(row_var.unit_factor), ['Der ''unit_factor'' der Variable ',row_var.var_name{1}, ' ist keine interpretierbare Zahl'])
+
+ if strcmp(run_var{1},'Z_el') % da die Impedanz eine Frequenz als Input benötigt, um berechnet zu werden
+ assert(~isnan(options.frequency),'To be able to calculate the impedance of the bearing, a frequency must be specified in Hz! \n e.g. '' frequency'', 1e4')
+ var_raw = eval(row_var.var{1});
+ plot_var_raw = abs(var_raw(options.frequency)) *row_var.unit_factor;
+ else
+ plot_var_raw = eval(row_var.var{1})*row_var.unit_factor; % hier wird der unit_factor mit berücksichtigt
+ end
+
+ if size(plot_var_raw,3)>1 % Die Graphen werden immer auf den ersten Wälzkörper bezogen geplottet
+ plot_var_raw = squeeze(plot_var_raw(:,:,1));
+ end
+
+
+ if isnan(obj.psi_calc)
+
+ [plot_psi,ind] = check_psi(obj.psi,options.plotStyle);
+ plot_var = plot_var_raw(:,ind);
+
+ else
+
+ % da die Funktion calculate keine Werte doppelt berechnet, wird hier der Vetor wieder vervollständigt
+ if strcmp(row_var.var_name,'delta_sx')
+ [plot_psi,ind] = check_psi(obj.psi_all,options.plotStyle);
+ plot_var = plot_var_raw(:,ind);
+ % da delta_sx die einzige Variable ist, die punktsymmetrisch gekürzt wird
+
+ plot_var((obj.psi_calc(ind)>pi) & (plot_psi>0)) = -plot_var((obj.psi_calc(ind)>pi) & (plot_psi>0));
+ plot_var((obj.psi_calc(ind)<pi) & (plot_psi<0)) = -plot_var((obj.psi_calc(ind)<pi) & (plot_psi<0));
+ elseif any(strcmp(row_var.var_name,{'C_el_single','R_el_single','C_el','R_el','Z_el'})) && obj.L.allBallsSameMaterial
+ [plot_psi,ind] = check_psi(obj.psi,options.plotStyle);
+ plot_var = plot_var_raw(:,ind);
+ else
+ [plot_psi,ind] = check_psi(obj.psi_all,options.plotStyle);
+ plot_var = plot_var_raw(:,ind);
+ end
+
+ end
+
+ switch options.plotStyle
+ case 'linear'
+ plotLinear(plot_var,plot_psi, row_var.title_var, {row_var.full_name{1};strcat(row_var.var_label{1}," in ",row_var.unit{1})})
+ case 'polar'
+ plotPolar(plot_var,plot_psi, row_var.title_var, row_var.unit{1})
+ end
+ end
+
+function plotLinear(vek,psi, pltTitle, pltYlabel)
+% plots a linear graph with all necessary adjustments
+ vek(isnan(vek)) = 0;
+ plot(psi,vek,'linewidth',1) %length(a.B.delta_s(1,:)
+ title(pltTitle)
+ xlabel("Cage angle in °")
+ ylabel(pltYlabel)
+ set(gca,'XTick',(-pi:2*pi/12:pi));
+ set(gca,'XTickLabel',(-180:30:181));
+ xlim([-pi pi])
+ hold off
+ grid on
+ if size(vek,1)==2
+ legend('inner contact','outer contact') %,'Position',[0.76 0.83 0.05 0.07]
+ end
+end
+
+function plotPolar(vek,ind, pltTitle, unit)
+% plots a polar graph with all necessary adjustments
+ vek_plot = [vek,vek(:,1)];
+ vek_plot(isnan(vek_plot)) = 0;
+ vek_max = max(vek_plot(1,:));
+
+ if vek_max == inf
+ temp_vek = vek_plot(~isinf(vek_plot) );
+ vek_max = max(temp_vek);
+ end
+
+ vek_min = min(vek_plot(1,:));
+
+ digitsRound = floor(log10(abs(vek_max)))+2;
+
+ if abs(digitsRound) == inf
+ vek_ax_max = vek_max;
+ vek_ax_min = vek_min;
+ else
+ vek_ax_max = round(vek_max,digitsRound,'significant');
+ vek_ax_min = round(vek_min,digitsRound,'significant');
+ end
+
+
+ if vek_ax_min>vek_min
+ vek_diff=round(vek_ax_max-vek_ax_min);
+ vek_ax_min = vek_ax_min-vek_diff/5;
+ end
+
+ if (vek_ax_max-vek_ax_min==0) && (vek_max~=0)
+ vek_diff = vek_max*1.1;
+ elseif (vek_ax_max-vek_ax_min==0) && (vek_max==0)
+ vek_diff = vek_max+0.1;
+ else
+ vek_diff = (vek_ax_max-vek_ax_min);
+ end
+
+
+ psi = [ind,ind(1)];
+
+ string1 = num2str(linspace(0,360-360/obj.L.Z,obj.L.Z)');
+ string2 = nan(obj.L.Z,2);
+ for run_RE = 1:obj.L.Z
+ string2(run_RE,:) = ' °';
+ end
+ thetaTickLabel = [string1,string2];
+
+ polarplot(psi,vek_plot);
+ ax = gca;
+ ax.ThetaDir = 'clockwise';
+ ax.RAxis.Color = 'black';
+ ax.ThetaZeroLocation = 'bottom';
+ ax.RLim = [vek_ax_min-vek_diff*5/6 vek_max+vek_diff/4];
+ ax.RTick = [vek_ax_min, vek_ax_min+vek_diff/4, vek_ax_min+vek_diff/2, vek_ax_min+vek_diff*3/4, vek_ax_min+vek_diff];
+ ax.RTickLabel = {strcat(num2str(vek_ax_min)," ",unit), strcat(num2str(vek_ax_min+vek_diff/4)," ",unit), strcat(num2str(vek_ax_min+vek_diff/2)," ",unit), strcat(num2str(vek_ax_min+vek_diff*3/4)," ",unit), strcat(num2str(vek_ax_max)," ",unit)};
+ ax.ThetaTick = linspace(0,360-360/obj.L.Z,obj.L.Z);
+ ax.ThetaTickLabel = {thetaTickLabel};
+ ax.RAxisLocation = 180;
+
+ title(pltTitle)
+ hold on
+ polarplot(psi,ones(1,length(psi))*vek_ax_min,'black');
+ if size(vek,1)==2
+ legend('inner contact','outer contact','') %,'Position',[0.83 0.07 0.05 0.07]
+ end
+ hold off
+end
+
+function [psi_checked,ind] = check_psi(psi,plotStyle)
+
+ psi(psi<0) = 2*pi + psi(psi<0);
+ psi_raw = psi;
+ psi = sort(psi);
+
+ ind = nan(1,length(psi));
+ for pos = 1: length(psi)
+ [~,~,ind(pos)] = intersect(psi(pos),psi_raw);
+ end
+
+ ind = [ind(psi>=0),ind(psi<0)];
+ if ~any(isnan(obj.psi_calc)) && ~(any(strcmp(row_var.var_name,{'C_el_single','R_el_single','C_el','R_el','Z_el'})) && obj.L.allBallsSameMaterial)
+ ind = obj.ind_psi_all(ind);
+ end
+
+ switch plotStyle
+ case 'linear'
+ psi_checked = [psi(psi>=pi)-2*pi,psi(psi<pi)];
+ ind = [ind(psi>=pi),ind(psi<pi)];
+ case 'polar'
+ psi_checked = psi;
+ end
+end
end
\ No newline at end of file
diff --git a/@BearImp/set.m b/@BearImp/set.m
index 474b8684afa74c36c3dfe76a2d250419351806ef..e512d212daf3ab4c26c56055da0ccf0d8c3aac16 100644
--- a/@BearImp/set.m
+++ b/@BearImp/set.m
@@ -1,36 +1,36 @@
-function set(obj, name, value)
-%Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
-%auch mit BearImp-Arrays bzw. -Matrizen
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- arguments
- obj BearImp
- name char
- value
- end
-
- assert(isscalar(value) || all(size(obj) == size(value)),...
- 'Entweder muss value skalar sein oder obj und value müssen die gleiche Dimension aufweisen')
-
- if isscalar(obj) && isscalar(value)
- obj.(name) = value;
- elseif isscalar(value)
- setObjArrayPropToSingleValue(obj, name, value)
- else
- setObjArrayPropToValueArray(obj, name, value)
- end
-end
-
-function setObjArrayPropToSingleValue(objArray, name, value)
- for ii = 1:numel(objArray)
- objArray(ii).(name) = value;
- end
-end
-
-function setObjArrayPropToValueArray(objArray, name, valueArray)
- for ii = 1:numel(objArray)
- objArray(ii).(name) = valueArray(ii);
- end
+function set(obj, name, value)
+%Weist einem Attribut eines BearImp-Objekts einen Wert zu. Funktioniert
+%auch mit BearImp-Arrays bzw. -Matrizen
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+ arguments
+ obj BearImp
+ name char
+ value
+ end
+
+ assert(isscalar(value) || all(size(obj) == size(value)),...
+ 'Entweder muss value skalar sein oder obj und value müssen die gleiche Dimension aufweisen')
+
+ if isscalar(obj) && isscalar(value)
+ obj.(name) = value;
+ elseif isscalar(value)
+ setObjArrayPropToSingleValue(obj, name, value)
+ else
+ setObjArrayPropToValueArray(obj, name, value)
+ end
+end
+
+function setObjArrayPropToSingleValue(objArray, name, value)
+ for ii = 1:numel(objArray)
+ objArray(ii).(name) = value;
+ end
+end
+
+function setObjArrayPropToValueArray(objArray, name, valueArray)
+ for ii = 1:numel(objArray)
+ objArray(ii).(name) = valueArray(ii);
+ end
end
\ No newline at end of file
diff --git a/@BearImp/setFitting.m b/@BearImp/setFitting.m
index 8085b66eb3efa15d661b7a0bc5f7169b857c87b0..63ff0f18b8b2bc30cc69ca82cc1264e4fb7cbf37 100644
--- a/@BearImp/setFitting.m
+++ b/@BearImp/setFitting.m
@@ -1,74 +1,74 @@
-function setFitting(obj,shaft,housing)
-%setFitting -> öffnet Auswahlmenüs, um Welle und Gehäuse auszuwählen
-%setFitting('drawing') -> lädt die in der Zeichnung angegebenen Werte
-%setFitting(shaft,housing) -> lädt die angegebene Welle & Gehäuse
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- arguments
- obj
- shaft char = ''
- housing char = ''
- end
- if numel(obj) > 1
- for ii = 1:numel(obj) % TODO: not working with matrices
- obj(ii).setFitting(shaft,housing);
- end
- return
- end
- if nargin == 2 && strcmp(shaft,'drawing')
- housing = 'drawing';
- end
-
- fittingTable = readtable(obj.inputDataTableFile,'Sheet','Fittings');
- if isempty(shaft) % keine Welle gegeben
- [index, ok] = listdlg('ListString',fittingTable.Shaft,'SelectionMode','Single','Name','Auswahl Welle','PromptString','Wählen Sie eine Welle: ');
- if ~ok, return, end
- setShaftTolerances(obj,fittingTable,index)
- else % Welle gegeben
- if strcmp(shaft,'drawing')
- shaft = [num2str(obj.L.d*1e3) ' drawing'];
- end
- index = strcmp(shaft,fittingTable.Shaft);
- if ~any(index)
- error('Wellenbezeichnung nicht gefunden')
- else
- setShaftTolerances(obj,fittingTable,index)
- end
- end
- assert(fittingTable(index,:).d == obj.L.d,'Welle passt nicht zum Lager')
- obj.L.shaft = obj.materials.(fittingTable(index,:).material_shaft{:});
- % TODO: Redundanz beseitigen ########
- if isempty(housing) % kein Gehäuse gegeben
- [index, ok] = listdlg('ListString',fittingTable.Housing,'SelectionMode','Single','Name','Auswahl Gehäuse','PromptString','Wählen Sie ein Gehäuse: ');
- if ~ok, return, end
- setHousingTolerances(obj,fittingTable,index)
- else % Welle gegeben
- if strcmp(housing,'drawing')
- housing = [num2str(obj.L.D*1e3) ' drawing'];
- end
- index = strcmp(housing,fittingTable.Housing);
- if ~any(index)
- error('Gehäusebezeichnung nicht gefunden')
- else
- setHousingTolerances(obj,fittingTable,index)
- end
- end
- assert(fittingTable(index,:).D == obj.L.D,'Gehäuse passt nicht zum Lager')
- obj.L.housing = obj.materials.(fittingTable(index,:).material_housing{:});
- obj.L = orderfields(obj.L);
-end
-
-function setShaftTolerances(obj,fittingTable,index)
- obj.L.ei_d = fittingTable(index,:).Delta_lower_shaft;
- obj.L.es_d = fittingTable(index,:).Delta_upper_shaft;
- obj.L.d_i = fittingTable(index,:).d_i;
-end
-
-function setHousingTolerances(obj,fittingTable,index)
- obj.L.EI_D = fittingTable(index,:).Delta_lower_housing;
- obj.L.ES_D = fittingTable(index,:).Delta_upper_housing;
- obj.L.D_G = fittingTable(index,:).D_G;
-end
+function setFitting(obj,shaft,housing)
+%setFitting -> öffnet Auswahlmenüs, um Welle und Gehäuse auszuwählen
+%setFitting('drawing') -> lädt die in der Zeichnung angegebenen Werte
+%setFitting(shaft,housing) -> lädt die angegebene Welle & Gehäuse
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+ arguments
+ obj
+ shaft char = ''
+ housing char = ''
+ end
+ if numel(obj) > 1
+ for ii = 1:numel(obj) % TODO: not working with matrices
+ obj(ii).setFitting(shaft,housing);
+ end
+ return
+ end
+ if nargin == 2 && strcmp(shaft,'drawing')
+ housing = 'drawing';
+ end
+
+ fittingTable = readtable(obj.inputDataTableFile,'Sheet','Fittings');
+ if isempty(shaft) % keine Welle gegeben
+ [index, ok] = listdlg('ListString',fittingTable.Shaft,'SelectionMode','Single','Name','Auswahl Welle','PromptString','Wählen Sie eine Welle: ');
+ if ~ok, return, end
+ setShaftTolerances(obj,fittingTable,index)
+ else % Welle gegeben
+ if strcmp(shaft,'drawing')
+ shaft = [num2str(obj.L.d*1e3) ' drawing'];
+ end
+ index = strcmp(shaft,fittingTable.Shaft);
+ if ~any(index)
+ error('Wellenbezeichnung nicht gefunden')
+ else
+ setShaftTolerances(obj,fittingTable,index)
+ end
+ end
+ assert(fittingTable(index,:).d == obj.L.d,'Welle passt nicht zum Lager')
+ obj.L.shaft = obj.materials.(fittingTable(index,:).material_shaft{:});
+ % TODO: Redundanz beseitigen ########
+ if isempty(housing) % kein Gehäuse gegeben
+ [index, ok] = listdlg('ListString',fittingTable.Housing,'SelectionMode','Single','Name','Auswahl Gehäuse','PromptString','Wählen Sie ein Gehäuse: ');
+ if ~ok, return, end
+ setHousingTolerances(obj,fittingTable,index)
+ else % Welle gegeben
+ if strcmp(housing,'drawing')
+ housing = [num2str(obj.L.D*1e3) ' drawing'];
+ end
+ index = strcmp(housing,fittingTable.Housing);
+ if ~any(index)
+ error('Gehäusebezeichnung nicht gefunden')
+ else
+ setHousingTolerances(obj,fittingTable,index)
+ end
+ end
+ assert(fittingTable(index,:).D == obj.L.D,'Gehäuse passt nicht zum Lager')
+ obj.L.housing = obj.materials.(fittingTable(index,:).material_housing{:});
+ obj.L = orderfields(obj.L);
+end
+
+function setShaftTolerances(obj,fittingTable,index)
+ obj.L.ei_d = fittingTable(index,:).Delta_lower_shaft;
+ obj.L.es_d = fittingTable(index,:).Delta_upper_shaft;
+ obj.L.d_i = fittingTable(index,:).d_i;
+end
+
+function setHousingTolerances(obj,fittingTable,index)
+ obj.L.EI_D = fittingTable(index,:).Delta_lower_housing;
+ obj.L.ES_D = fittingTable(index,:).Delta_upper_housing;
+ obj.L.D_G = fittingTable(index,:).D_G;
+end
\ No newline at end of file
diff --git a/@BearImp/setLube.m b/@BearImp/setLube.m
index c3daf9d5ea0b459c864e45d99ee81e671e0e8906..0398f406c37a0f62f9fb2f3f77013a71b7e55fa4 100644
--- a/@BearImp/setLube.m
+++ b/@BearImp/setLube.m
@@ -1,34 +1,34 @@
-function setLube(obj,name)
-%setLube -> öffnet ein Auswahlmenü, um einen Schmierstoff auszuwählen
-%setLube(name) -> lädt den angegebenen Schmierstoff
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler
-
- arguments
- obj
- name char = ''
- end
-
- if numel(obj) > 1
- for ii = 1:numel(obj) % TODO: not working with matrices
- obj(ii).setLube(name);
- end
- return
- end
-
- lubeTable = readtable(obj.inputDataTableFile,'Sheet','Lubricants');
- if isempty(name) % kein Schmierstoff gegeben
- [selection, ok] = listdlg('ListString',lubeTable.Labelling,'SelectionMode','Single','Name','Schmierstoffauswahl','PromptString','Wählen Sie einen Schmierstoff: ');
- if ~ok, return, end
- S = table2struct(lubeTable(selection,:));
- else % Schmierstoff gegeben
- index = strcmp(name,lubeTable.Labelling);
- if ~any(index)
- error('Schmierstoffname nicht gefunden')
- else
- S = table2struct(lubeTable(index,:));
- end
- end
- obj.S = S;
+function setLube(obj,name)
+%setLube -> öffnet ein Auswahlmenü, um einen Schmierstoff auszuwählen
+%setLube(name) -> lädt den angegebenen Schmierstoff
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler
+
+ arguments
+ obj
+ name char = ''
+ end
+
+ if numel(obj) > 1
+ for ii = 1:numel(obj) % TODO: not working with matrices
+ obj(ii).setLube(name);
+ end
+ return
+ end
+
+ lubeTable = readtable(obj.inputDataTableFile,'Sheet','Lubricants');
+ if isempty(name) % kein Schmierstoff gegeben
+ [selection, ok] = listdlg('ListString',lubeTable.Labelling,'SelectionMode','Single','Name','Schmierstoffauswahl','PromptString','Wählen Sie einen Schmierstoff: ');
+ if ~ok, return, end
+ S = table2struct(lubeTable(selection,:));
+ else % Schmierstoff gegeben
+ index = strcmp(name,lubeTable.Labelling);
+ if ~any(index)
+ error('Schmierstoffname nicht gefunden')
+ else
+ S = table2struct(lubeTable(index,:));
+ end
+ end
+ obj.S = S;
end
\ No newline at end of file
diff --git a/Compare2OnlineCalc.mlx b/Compare2OnlineCalc.mlx
index 6089e64924b3f8d71755c46fbc52d27be5f65b72..7087a512d984830ef55552fb366aa8da1151aef7 100644
Binary files a/Compare2OnlineCalc.mlx and b/Compare2OnlineCalc.mlx differ
diff --git a/PlotLabel.csv b/PlotLabel.csv
index c3d0f707b16ebea3de1ef2af7747d9da1bc59f59..a89755a70394f0be7155c29ae7099ac1b86fa001 100644
--- a/PlotLabel.csv
+++ b/PlotLabel.csv
@@ -1,15 +1,17 @@
-var_name;var_label;var;full_name;unit;unit_factor;title_var
-delta_sy;\delta_{sy};obj.B.delta_s(1,:);vertial shaft displacement;\mum;1.00E+06;Shaft displacement in vertical direction over cage angle
-delta_sx;\delta_{sx};obj.B.delta_s(2,:);horizontal shaft displacement;\mum;1.00E+06;Shaft displacement in horizontal direction over cage angle
-Q;Q;obj.B.Q;load of single RE;kN;1.00E-03;Load of single RE over cage angle
-s;s;obj.B.s;intersection of ball and raceway;\mum;1.00E+06;Intersection of ball and raceway over cage angle
-h_0;h_0;obj.H.h_0;central film thickness;\mum;1.00E+06;Central film thickness over cage angle
-C_el_single;C_{el,single};obj.C.C_el_single;capacity of sinlge RE;pF;1.00E+12;Capacity of single RE over cage angle
-R_el_single;R_{el,single};obj.C.R_el_single;resistance of sinlge RE;k\Omega;1.00E-03;Resistance of single RE over cage angle
-C_el;C_{el};obj.C.C_el;total capacity of bearing;pF;1.00E+12;Total capacity of bearing over cage angle
-R_el;R_{el};obj.C.R_el;total resistance of bearing;k\Omega;1.00E-03;Total resistance of bearing over cage angle
-Z_el;|Z_{el}|;obj.Z.Z_el;total absolute impedance of bearing;k\Omega;1.00E-03;Total impedance of bearing over cage angle
-p;p;obj.H.p;contact pressure;Gpa;1.00E-09;Contact pressure over cage angle
-C_Hertz;C_{Hertz};obj.C.C_Hertz;capacity of Hertzian contact;pF;1.00E+12;Capacity of Hertzian contact over cage angle
-C_out;C_{out};obj.C.C_out;capacity of outer RE area;pF;1.00E+12;Capacity of outer RE area over cage angle
-intersection_RE;\delta_{inters};obj.B.intersectionBall;intersection of RE and raceway;\mum;1.00E+06;Intersection of ball and raceway over cage angle without film
+var_name;var_label;var;full_name;unit;unit_factor;title_var
+delta_sy;\delta_{sy};obj.B.delta_s(2,:);Vertial shaft displacement;\mum;1.00E+06;Shaft displacement in vertical direction over cage angle
+delta_sx;\delta_{sx};obj.B.delta_s(1,:);Horizontal shaft displacement;\mum;1.00E+06;Shaft displacement in horizontal direction over cage angle
+delta_sz;\delta_{sz};obj.B.delta_s(3,:);Axial shaft displacement;\mum;1.00E+06;Shaft displacement in axial direction over cage angle
+Q;Q;obj.B.Q;Load of single RE;kN;1.00E-03;Load of single RE over cage angle
+s;s;obj.B.s;Gap between ball and raceway;\mum;1.00E+06;Gap between ball and raceway over cage angle
+h_0;h_0;obj.H.h_0;Central film thickness;\mum;1.00E+06;Central film thickness over cage angle
+C_el_single;C_{el,single};obj.C.C_el_single;Capacitance of sinlge RE;pF;1.00E+12;Capacitance of single RE over cage angle
+R_el_single;R_{el,single};obj.C.R_el_single;Resistance of sinlge RE;k\Omega;1.00E-03;Resistance of single RE over cage angle
+C_el;C_{el};obj.C.C_el;Total capacitance of bearing;pF;1.00E+12;Total capacitance of bearing over cage angle
+R_el;R_{el};obj.C.R_el;Total resistance of bearing;k\Omega;1.00E-03;Total resistance of bearing over cage angle
+Z_el;|Z_{el}|;obj.Z.Z_el;Total absolute impedance of bearing;k\Omega;1.00E-03;Total impedance of bearing over cage angle
+p;p;obj.H.p;Contact pressure;Gpa;1.00E-09;Contact pressure over cage angle
+C_Hertz;C_{Hertz};obj.C.C_Hertz;Capacitance of Hertzian contact;pF;1.00E+12;Capacitance of Hertzian contact over cage angle
+C_out;C_{out};obj.C.C_out;Capacitance of outer RE area;pF;1.00E+12;Capacitance of outer RE area over cage angle
+intersection_RE;\delta_{inters};obj.B.intersectionBall(:,:,1);Intersection of RE and raceway;\mum;1.00E+06;Intersection of ball and raceway over cage angle without film
+alpha;\alpha;obj.B.alpha;Load angle;�;57;Load angle over cage angle
diff --git a/possibleMethods.m b/possibleMethods.m
index 18d4f07137efa63cf0804bb09288ffd25650c779..e34c17c0d200706d4beed3d41273e967f937cfe6 100644
--- a/possibleMethods.m
+++ b/possibleMethods.m
@@ -1,112 +1,112 @@
-classdef possibleMethods
-% Gibt Auskunft über mögliche Berechnungsmethoden der einzelnen
-% Berechnungsabschnitte (T, R, G, B, H, C) und enthält eine Funktion zum
-% überprüfen des Methodenstructs
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Steffen Puchtler, Julius van der Kuip
-
- methods (Static) % Diese Klasse enthält ausschließlich statische Methoden, sodass kein Objekt initialisiert werden muss
- function s = T
- s = {'linear','Vogel'};
- end
- function s = R
- s = {};
- end
- function s = G
- s = {};
- end
- function s = B
- s = {'static','dynamic'};
- end
- function s = H
- s = {'Hamrock/Dowson','Moes'};
- end
- function s = C
- s.unloadedRE = { 'neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
- s.outsideArea = {'k-factor','neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
- end
-
- function s = Default
- % Gibt ein Rechenmethoden-Struct mit den Default-Methoden aus
- s.T = 'Vogel';
- s.B = 'static';
- s.H = 'Hamrock/Dowson';
- s.C.unloadedRE = 'semianalytisch3D';
- s.C.outsideArea = 'semianalytisch3D';
- end
-
- function method = addDefault(method)
- % Ergänzt ein unvollständiges Rechenmethoden-Struct mit den
- % Default-Werten
- arguments
- method (1,1) struct
- end
- defaultMethod = possibleMethods.Default;
- for ii = fields(defaultMethod)'
- if ~isfield(method,ii{1}) % Falls Feld ii nicht besetzt,
- method.(ii{1}) = defaultMethod.(ii{1}); % überschreibe mit default
- elseif isstruct(defaultMethod.(ii{1})) % falls zweite Ebene existiert
- for jj = fields(defaultMethod.(ii{1}))'
- if ~isfield(method.(ii{1}),jj{1}) % Falls Feld ii in jj nicht besetzt
- method.(ii{1}).(jj{1}) = defaultMethod.(ii{1}).(jj{1});
- end
- end
- end
- end
- end
-
- function success = check(method,throwError)
- % prüft, ob das Struct method ausschließlich valide Feld- und Methodennamen enthält
- % optional: throwError = true gibt error-Meldungen aus
- arguments
- method (1,1) struct
- throwError (1,1) logical = false
- end
- success = true;
- for ii = fields(method)'
- if ~any(strcmp(ii,methods(possibleMethods)')) % Prüft Feldname in erster Ebene
- success = false;
- if throwError
- error(['%s ist kein gültiger Feldname im Methoden-Struct\n',...
- 'Möglich an dieser Stelle: %s'],...
- ii{1}, strjoin(methods(possibleMethods)') )
- end
- break
- end
- if isstruct(method.(ii{1})) % falls zweite Ebene existiert
- for jj = fields(method.(ii{1}))'
- if ~any(strcmp(jj,fields(possibleMethods.(ii{1}))')) % Prüft Feldname in zweiter Ebene
- success = false;
- if throwError
- error(['%s is kein gültiger Feldname für den Berechnungsabschnitt %s\n',...
- 'Möglich an dieser Stelle: %s'],...
- jj{1},ii{1},strjoin(fields(possibleMethods.(ii{1}))') )
- end
- break
- end
- if ~any(strcmp(method.(ii{1}).(jj{1}),possibleMethods.(ii{1}).(jj{1}))) % Prüft Methodenname in zweiter Ebene
- success = false;
- if throwError
- error(['%s ist kein gültiger Methodenname zur Berechnung von %s im Berechnungsabschnitt %s\n',...
- 'Möglich an dieser Stelle: %s'],...
- method.(ii{1}).(jj{1}), jj{1}, ii{1}, strjoin(possibleMethods.(ii{1}).(jj{1})) )
- end
- break
- end
- end
- else
- if ~any(strcmp(method.(ii{1}),possibleMethods.(ii{1}))) % Prüft Methodenname in erster Ebene
- success = false;
- if throwError
- error(['%s ist kein gültiger Methodenname für den Berechnungsabschnitt %s\n',...
- 'Möglich an dieser Stelle: %s'],...
- method.(ii{1}), ii{1}, strjoin(possibleMethods.(ii{1})) )
- end
- break
- end
- end
- end
- end
- end
+classdef possibleMethods
+% Gibt Auskunft über mögliche Berechnungsmethoden der einzelnen
+% Berechnungsabschnitte (T, R, G, B, H, C) und enthält eine Funktion zum
+% überprüfen des Methodenstructs
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Steffen Puchtler, Julius van der Kuip
+
+ methods (Static) % Diese Klasse enthält ausschließlich statische Methoden, sodass kein Objekt initialisiert werden muss
+ function s = T
+ s = {'linear','Vogel'};
+ end
+ function s = R
+ s = {};
+ end
+ function s = G
+ s = {};
+ end
+ function s = B
+ s = {'static','dynamic'};
+ end
+ function s = H
+ s = {'Hamrock/Dowson','Moes'};
+ end
+ function s = C
+ s.unloadedRE = { 'neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
+ s.outsideArea = {'k-factor','neglect','Leander_Parallel','Leander_Radial','LeanderSteffen','TobiasSteffen_Kugelfläche','TobiasSteffen_Laufbahnfläche','semianalytisch3D'};
+ end
+
+ function s = Default
+ % Gibt ein Rechenmethoden-Struct mit den Default-Methoden aus
+ s.T = 'Vogel';
+ s.B = 'static';
+ s.H = 'Hamrock/Dowson';
+ s.C.unloadedRE = 'semianalytisch3D';
+ s.C.outsideArea = 'semianalytisch3D';
+ end
+
+ function method = addDefault(method)
+ % Ergänzt ein unvollständiges Rechenmethoden-Struct mit den
+ % Default-Werten
+ arguments
+ method (1,1) struct
+ end
+ defaultMethod = possibleMethods.Default;
+ for ii = fields(defaultMethod)'
+ if ~isfield(method,ii{1}) % Falls Feld ii nicht besetzt,
+ method.(ii{1}) = defaultMethod.(ii{1}); % überschreibe mit default
+ elseif isstruct(defaultMethod.(ii{1})) % falls zweite Ebene existiert
+ for jj = fields(defaultMethod.(ii{1}))'
+ if ~isfield(method.(ii{1}),jj{1}) % Falls Feld ii in jj nicht besetzt
+ method.(ii{1}).(jj{1}) = defaultMethod.(ii{1}).(jj{1});
+ end
+ end
+ end
+ end
+ end
+
+ function success = check(method,throwError)
+ % prüft, ob das Struct method ausschließlich valide Feld- und Methodennamen enthält
+ % optional: throwError = true gibt error-Meldungen aus
+ arguments
+ method (1,1) struct
+ throwError (1,1) logical = false
+ end
+ success = true;
+ for ii = fields(method)'
+ if ~any(strcmp(ii,methods(possibleMethods)')) % Prüft Feldname in erster Ebene
+ success = false;
+ if throwError
+ error(['%s ist kein gültiger Feldname im Methoden-Struct\n',...
+ 'Möglich an dieser Stelle: %s'],...
+ ii{1}, strjoin(methods(possibleMethods)') )
+ end
+ break
+ end
+ if isstruct(method.(ii{1})) % falls zweite Ebene existiert
+ for jj = fields(method.(ii{1}))'
+ if ~any(strcmp(jj,fields(possibleMethods.(ii{1}))')) % Prüft Feldname in zweiter Ebene
+ success = false;
+ if throwError
+ error(['%s is kein gültiger Feldname für den Berechnungsabschnitt %s\n',...
+ 'Möglich an dieser Stelle: %s'],...
+ jj{1},ii{1},strjoin(fields(possibleMethods.(ii{1}))') )
+ end
+ break
+ end
+ if ~any(strcmp(method.(ii{1}).(jj{1}),possibleMethods.(ii{1}).(jj{1}))) % Prüft Methodenname in zweiter Ebene
+ success = false;
+ if throwError
+ error(['%s ist kein gültiger Methodenname zur Berechnung von %s im Berechnungsabschnitt %s\n',...
+ 'Möglich an dieser Stelle: %s'],...
+ method.(ii{1}).(jj{1}), jj{1}, ii{1}, strjoin(possibleMethods.(ii{1}).(jj{1})) )
+ end
+ break
+ end
+ end
+ else
+ if ~any(strcmp(method.(ii{1}),possibleMethods.(ii{1}))) % Prüft Methodenname in erster Ebene
+ success = false;
+ if throwError
+ error(['%s ist kein gültiger Methodenname für den Berechnungsabschnitt %s\n',...
+ 'Möglich an dieser Stelle: %s'],...
+ method.(ii{1}), ii{1}, strjoin(possibleMethods.(ii{1})) )
+ end
+ break
+ end
+ end
+ end
+ end
+ end
end
\ No newline at end of file
diff --git a/testProgram.m b/testProgram.m
index 844639699a76487e723e8b831a4a0cefe0cab72a..c72da0db77611958690d1c4c597829b3a1aa425f 100644
--- a/testProgram.m
+++ b/testProgram.m
@@ -1,144 +1,144 @@
-% Skript, um die Funktionsfähigkeit von BearImp zu testen. Muss vor jedem
-% Merge ausgeführt werden. Neue Berechnungsmethoden etc. müssen hier
-% ergänzt werden.
-
-% pmd Berechnungstool Lagerimpedanz
-% Autor: Julius van der Kuip
-
-clear all
-close all
-%% 1. Testlauf 'default' (a) [Standartfall Stahllager mit Voreinstellungen]
-
-% Berechnung
-a = BearImp('default');
-a.calculate
-
-a.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
- 'figureName',['Test1_Stahl_default: ' a.L.Labelling ' bei F_r=' num2str(a.F_r),...
- ' N und omega= ' num2str(a.omega), ' 1/s und T_Oil= ' num2str(a.T_Oil) ' °C'])
-
-msgText = '1. Testlauf erfolgreich abgeschlosen';
-msg = msgbox(msgText,'Programm Update');
-
-%% 2. Testlauf (b) [werden ohne obj.calculate trotzdem bei korrektem psi die richtigen Werte berechnet?']
-
-% Berechnung
-b = BearImp('default');
-
-b.calcLub
-b.calcClear
-b.calcGeo
-b.calcLoad
-b.calcFilm
-b.calcCap
-b.calcImp
-
-b.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
- 'figureName',['Test2_Stahl_ohneCalculate: ' b.L.Labelling ' bei F_r=' num2str(b.F_r), ...
- ' N und omega= ' num2str(b.omega), ' 1/s und T_Oil= ' num2str(b.T_Oil) ' °C'])
-
-msgText = {msgText;''; '2. Testlauf erfolgreich abgeschlosen!';'Vergleichen Sie die Verläufe des 1. Testlaufs mit den Verläufen des 2. Testlaufs.';'-> Diese sollten identisch aussehen!'};
-msgbox(msgText,'Programm Update','replace' );
-
-%% 3. Testlauf (c) [Mixlager mit default-Einstellungen]
-% Berechnung
-c = BearImp('default');
-c.setBearing("6205-C3 mix")
-c.calculate
-
-assert(length(c.psi) == length(c.psi_calc), 'psi_calc entspricht nicht dem optimalen Winkelvektor! Es ist die Funktion ''find_nessecary_psi'' in caculate zu überprüfen.')
-
-c.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
- 'figureName',['Test3_Mix_default: ' c.L.Labelling ' bei F_r=' num2str(c.F_r), ...
- ' N und omega= ' num2str(c.omega), ' 1/s und T_Oil= ' num2str(c.T_Oil) ' °C'])
-
-msgText{6} = '';
-msgText{7} = '3. Testlauf erfolgreich abgeschlosen!';
-msgbox(msgText,'Programm Update', 'replace');
-
-%% 4. Testlauf (d) [Mix-Lager mit wenigeren Abtastpunkten]
-% Berechnung
-d = BearImp;
-
-d.F_r = 3000; %Bereich zwischen 920-1020
-d.F_a = 0;
-d.omega=2000*2*pi/60;
-d.resolutionPerRotation=10;
-d.T_Oil = 50;
-d.setBearing("6205-C3 mix")
-d.setLube("FVA Referenz-Öl III")
-
-d.calculate
-
-d.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
- 'figureName',['Test4_Mix_wenigerPunkte: ' d.L.Labelling ' bei F_r=' num2str(d.F_r), ...
- ' N und omega= ' num2str(d.omega), ' 1/s und T_Oil= ' num2str(d.T_Oil) ' °C'])
-
-msgText{8} = '';
-msgText{9} = '4. Testlauf erfolgreich abgeschlosen!';
-msgbox(msgText,'Programm Update','replace');
-
-%% 5. Testlauf (e,f) [werden Fehlermeldungen richtig ausgeführt?]
-
-e = BearImp('default');
-e.psi = [0, 20, 30, 90, 130, 180, 201, 200, 270, 360]/360*2*pi;
-try
- e.calculate
- catch ME1
- error('Es wurden fälschlicher Weise nicht alle benötigten Werte für die Gesamtkapazitätsberechnung berechnet! (siehe ME1)')
-end
-try
-check_all_necessary_psi(e.psi,e.psi,e.L.Z,true)
-catch ME2
- strMsg2 = 'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!';
- assert(strcmp(strMsg2,ME2.message), 'Obwohl nicht alle Werte für Gesamtkapazität berechnet wurden, wird kein Error angezeigt! (siehe ME2)')
-end
-
-f = BearImp('default');
-f.calculate
-try
-check_all_necessary_psi(f.psi_calc,f.psi,f.L.Z,true)
-catch ME3
- error('psi_calc ermittelt nicht alle benötigten Winkel zur Berechnung der Gesamtkapazität! (siehe ME3)')
-end
-
-try
- f.F_r = 200;
- f.calcImp
-
-catch ME4
- strMsg4 = 'Kapazitäten noch nicht berechnet';
- assert(strcmp(strMsg4,ME4.message), 'up2date beinhaltet fehler, da nach Veränderung der radialen Kraft nicht alle benötigten Werte berechnet wurden und das Programm es trotzdem zulässt! (siehe ME4) ')
-end
-
-msgText{10} = '';
-msgText{11} = '5. Testlauf erfolgreich abgeschlosen!';
-msgText{12} = 'Fehlermeldungen werden richtig ausgeführt und schützen den Benutzer vor Fehlberechnungen.';
-msgText{13} = '';
-msgText{14} = 'Das TestProgram ist durchgeführt!';
-msgbox(msgText,'Programm Update', 'replace');
-
-%% Helper functions
-
-function check_all_necessary_psi(psi,obj_psi,Z,symCase)
-
- if symCase
- psi_calc = [psi, 2*pi - psi];
- else
- psi_calc = psi;
- end
-
- for run_C = 1 : length(obj_psi)
-
- psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj_psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
-
- psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch %*0.9999999
-
- psiNeedForIndCapa = round(psiNeedForIndCapa,8);
-
- [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-6);
-
- %Check, if all nessecary C_el_single are calculated
- assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!')
- end
+% Skript, um die Funktionsfähigkeit von BearImp zu testen. Muss vor jedem
+% Merge ausgeführt werden. Neue Berechnungsmethoden etc. müssen hier
+% ergänzt werden.
+
+% pmd Berechnungstool Lagerimpedanz
+% Autor: Julius van der Kuip
+
+clear all
+close all
+%% 1. Testlauf 'default' (a) [Standartfall Stahllager mit Voreinstellungen]
+
+% Berechnung
+a = BearImp('default');
+a.calculate
+
+a.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
+ 'figureName',['Test1_Stahl_default: ' a.L.Labelling ' bei F_r=' num2str(a.F_r),...
+ ' N und omega= ' num2str(a.omega), ' 1/s und T_Oil= ' num2str(a.T_Oil) ' °C'])
+
+msgText = '1. Testlauf erfolgreich abgeschlosen';
+msg = msgbox(msgText,'Programm Update');
+
+%% 2. Testlauf (b) [werden ohne obj.calculate trotzdem bei korrektem psi die richtigen Werte berechnet?']
+
+% Berechnung
+b = BearImp('default');
+
+b.calcLub
+b.calcClear
+b.calcGeo
+b.calcLoad
+b.calcFilm
+b.calcCap
+b.calcImp
+
+b.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
+ 'figureName',['Test2_Stahl_ohneCalculate: ' b.L.Labelling ' bei F_r=' num2str(b.F_r), ...
+ ' N und omega= ' num2str(b.omega), ' 1/s und T_Oil= ' num2str(b.T_Oil) ' °C'])
+
+msgText = {msgText;''; '2. Testlauf erfolgreich abgeschlosen!';'Vergleichen Sie die Verläufe des 1. Testlaufs mit den Verläufen des 2. Testlaufs.';'-> Diese sollten identisch aussehen!'};
+msgbox(msgText,'Programm Update','replace' );
+
+%% 3. Testlauf (c) [Mixlager mit default-Einstellungen]
+% Berechnung
+c = BearImp('default');
+c.setBearing("6205-C3 mix")
+c.calculate
+
+assert(length(c.psi) == length(c.psi_calc), 'psi_calc entspricht nicht dem optimalen Winkelvektor! Es ist die Funktion ''find_nessecary_psi'' in caculate zu überprüfen.')
+
+c.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
+ 'figureName',['Test3_Mix_default: ' c.L.Labelling ' bei F_r=' num2str(c.F_r), ...
+ ' N und omega= ' num2str(c.omega), ' 1/s und T_Oil= ' num2str(c.T_Oil) ' °C'])
+
+msgText{6} = '';
+msgText{7} = '3. Testlauf erfolgreich abgeschlosen!';
+msgbox(msgText,'Programm Update', 'replace');
+
+%% 4. Testlauf (d) [Mix-Lager mit wenigeren Abtastpunkten]
+% Berechnung
+d = BearImp;
+
+d.F_r = 3000; %Bereich zwischen 920-1020
+d.F_a = 0;
+d.omega=2000*2*pi/60;
+d.resolutionPerRotation=10;
+d.T_Oil = 50;
+d.setBearing("6205-C3 mix")
+d.setLube("FVA Referenz-Öl III")
+
+d.calculate
+
+d.plot({'delta_sy','delta_sx','Q','s','h_0','C_el_single','R_el_single','C_el','R_el','Z_el'},'frequency',30e3,...
+ 'figureName',['Test4_Mix_wenigerPunkte: ' d.L.Labelling ' bei F_r=' num2str(d.F_r), ...
+ ' N und omega= ' num2str(d.omega), ' 1/s und T_Oil= ' num2str(d.T_Oil) ' °C'])
+
+msgText{8} = '';
+msgText{9} = '4. Testlauf erfolgreich abgeschlosen!';
+msgbox(msgText,'Programm Update','replace');
+
+%% 5. Testlauf (e,f) [werden Fehlermeldungen richtig ausgeführt?]
+
+e = BearImp('default');
+e.psi = [0, 20, 30, 90, 130, 180, 201, 200, 270, 360]/360*2*pi;
+try
+ e.calculate
+ catch ME1
+ error('Es wurden fälschlicher Weise nicht alle benötigten Werte für die Gesamtkapazitätsberechnung berechnet! (siehe ME1)')
+end
+try
+check_all_necessary_psi(e.psi,e.psi,e.L.Z,true)
+catch ME2
+ strMsg2 = 'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!';
+ assert(strcmp(strMsg2,ME2.message), 'Obwohl nicht alle Werte für Gesamtkapazität berechnet wurden, wird kein Error angezeigt! (siehe ME2)')
+end
+
+f = BearImp('default');
+f.calculate
+try
+check_all_necessary_psi(f.psi_calc,f.psi,f.L.Z,true)
+catch ME3
+ error('psi_calc ermittelt nicht alle benötigten Winkel zur Berechnung der Gesamtkapazität! (siehe ME3)')
+end
+
+try
+ f.F_r = 200;
+ f.calcImp
+
+catch ME4
+ strMsg4 = 'Kapazitäten noch nicht berechnet';
+ assert(strcmp(strMsg4,ME4.message), 'up2date beinhaltet fehler, da nach Veränderung der radialen Kraft nicht alle benötigten Werte berechnet wurden und das Programm es trotzdem zulässt! (siehe ME4) ')
+end
+
+msgText{10} = '';
+msgText{11} = '5. Testlauf erfolgreich abgeschlosen!';
+msgText{12} = 'Fehlermeldungen werden richtig ausgeführt und schützen den Benutzer vor Fehlberechnungen.';
+msgText{13} = '';
+msgText{14} = 'Das TestProgram ist durchgeführt!';
+msgbox(msgText,'Programm Update', 'replace');
+
+%% Helper functions
+
+function check_all_necessary_psi(psi,obj_psi,Z,symCase)
+
+ if symCase
+ psi_calc = [psi, 2*pi - psi];
+ else
+ psi_calc = psi;
+ end
+
+ for run_C = 1 : length(obj_psi)
+
+ psiNeedForIndCapa = ((1:Z)' -1) / Z * 2*pi + obj_psi(run_C); % Bsp. Ergebnis: [0:2*pi/L.Z:2*pi] im Abstand der WK und dem gewünschten Offset von psi
+
+ psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) = psiNeedForIndCapa(psiNeedForIndCapa >= 2*pi - eps(3*pi)) - 2*pi; %alle Winkel kleiner 2*pi, da periodisch %*0.9999999
+
+ psiNeedForIndCapa = round(psiNeedForIndCapa,8);
+
+ [~, ind_single_psi] = ismembertol(psiNeedForIndCapa, psi_calc,1e-6);
+
+ %Check, if all nessecary C_el_single are calculated
+ assert(~any(ind_single_psi==0),'Es wurden nicht alle benötigten Einzelnkapazitäten berechnet, um die Gesamtkapazität für eine Position zu berechnen!')
+ end
end
\ No newline at end of file