基于MATLAB实现的实现设置滚都轴承的参数，然后计算滚都轴承的刚度矩阵，位移和受力+使用说明文档.zip

%clear all %clc %close all %% Initialise Global Varibles n = 3/2; k = 3.5e9; % load-deflection factor for hertzian theory nb = 9; % the number of rolling elements %rb = 4.7625e-3; % the radius of the rolling element rb = 3.97; A0 = 200e-6; % distance between the inner and outer curvature of the groove when unloaded alpha0 = 0/180*pi; %unloaded contact angle rcl = 0.5e-6; % radial clearance r_curve = 19.65e-3; % Radius of the inner raceway groove curvature centre R0 = rb + r_curve; % Radius of the outer raceway groove curvature centre Load = [0 -5000 0 0 0]'; %% Define the Defect phi_axis = linspace(0, 2*pi, 360/0.01); defect_pos = 180*pi/180; % the location of the centre of the defect defect_size = [0 15.8 55.8]*pi/180; % the length of the defect in radians defect_depth = 100e-6; % the depth of the defect phi_0 = 2*pi*[0:nb-1]/nb; % the initial loctions of the rolling elements phi_cage = 0:2*pi/(360*nb):2*pi/nb*(1-1/360); % the angle of rotation of the cage %% Initialise Load Matrices Q_force = zeros(length(phi_cage),nb,length(defect_size)); % The contact force on the j element delta = zeros(length(phi_cage),nb,length(defect_size)); % the radial displace for each rolling element delta_g = zeros(length(phi_cage),5,length(defect_size)); % the global diplacement and moment matrix Fxs = zeros(length(phi_cage),nb,length(defect_size)); % the force in the x direction acting on the jth element Fys = zeros(length(phi_cage),nb,length(defect_size)); % the force in the y direction acting on the jth element Fzs = zeros(length(phi_cage),nb,length(defect_size)); % the force in the x direction acting on the jth element Mxs = zeros(length(phi_cage),nb,length(defect_size)); % the moment about the x axis caused by the jth element Mys = zeros(length(phi_cage),nb,length(defect_size)); % the moment about the y axis caused by the jth element % Stiffness matrices % kxx is the stiffness coeff of x axis, kxy is the stiffness coeff of x % coupled with y, and so on. While kxrx is the stiffness coeff. of x % coupled with the rotate of the x axis, while krxrx is the stiffness coeff % of the rotation about the x axis, and so on. kxx = zeros(length(phi_cage),length(defect_size)); kxy = zeros(length(phi_cage),length(defect_size)); kxz = zeros(length(phi_cage),length(defect_size)); kxrx = zeros(length(phi_cage),length(defect_size)); kxry = zeros(length(phi_cage),length(defect_size)); kyy = zeros(length(phi_cage),length(defect_size)); kyz = zeros(length(phi_cage),length(defect_size)); kyrx = zeros(length(phi_cage),length(defect_size)); kyry = zeros(length(phi_cage),length(defect_size)); kzz = zeros(length(phi_cage),length(defect_size)); kzrx = zeros(length(phi_cage),length(defect_size)); kzry = zeros(length(phi_cage),length(defect_size)); krxrx = zeros(length(phi_cage),length(defect_size)); krxry = zeros(length(phi_cage),length(defect_size)); kryry = zeros(length(phi_cage),length(defect_size)); phi_m = zeros(length(phi_cage),nb); % the full list of element position throughout the simulation Cd = zeros(length(phi_axis),1); %% Initialize Progress Bar %progress_bar = waitbar(0,'Calculating Matrices'); %nruns = length(phi_cage)*length(defect_size); %number of runs equal the number of different defect sizes time by the number a cage movements run = 1; %% Calculate the displacments, Forces and Stiffness Matrices set(gcf,'Position',[1 29 1685 955]); % for each defect size calculate the displacements, forces and Stiffness % matrices for i= 1:length(defect_size) %plot the defect depth against position on the outer raceway Cd(:,i) = profile_square(phi_axis,defect_size(i),defect_pos,defect_depth,R0,rb); subplot(2,2,i); plot(phi_axis*180/pi,Cd(:,i)*1e6); grid on; ylim([0 60]); xlim([0 360]); set(gca,'XTick',[0:45:360]); ylabel('depth (\mum)'); xlabel('\phi(\circ)'); for j = 1:length(phi_cage) %waitbar(run/nruns,progress_bar); % the position of the rolling elements with in the bearing %{ if(phi_cage(j)< 10*pi/180 && phi_cage(j)< 30*pi/180) if(i>2) if(Load(1)*-1<550) k = 1e9; elseif(Load(1)*-1<1550) k = 2e9; else k=7e9; end else k=7e9; end end %} phi = phi_0 + phi_cage(j); phi_m(j,:) = phi; % defect depth at rolling element positions Cdi = interp1(phi_axis,Cd(:,i),phi); % initial guess for the displacment x0 = [0 -5e-6 0 0 0]'; % found the solution for the displacement of raceways delta_g(j,:,i) = FindBallBearingLoad(phi,Load,rcl,Cdi,r_curve,alpha0,A0,k,x0,n); x = delta_g(j,:,i); % Contact deformations delta_r = x(1)*cos(phi)+x(2)*sin(phi)-rcl-Cdi*cos(alpha0); delta_z = x(3)+r_curve*(x(4)*sin(phi)-x(5)*cos(phi))-Cdi*sin(alpha0); delta_rs = A0*cos(alpha0)+delta_r; delta_zs = A0*sin(alpha0)+delta_z; alpha = atan2(delta_zs,delta_rs); A = sqrt(delta_zs.^2+delta_rs.^2); delta(j,:,i) = max(A-A0,0); I = find(delta(j,:,i)>0); D_zs(j,:,i) = delta_zs; D_rs(j,:,i) = delta_rs; Alpha(j,:,i) = alpha; % Contact forces and moments Q_force(j,:,i) = k*delta(j,:,i).^n; %delta Fxs(j,:,i) = Q_force(j,:,i).*cos(alpha).*cos(phi); Fys(j,:,i) = Q_force(j,:,i).*cos(alpha).*sin(phi); Fzs(j,:,i) = Q_force(j,:,i).*sin(alpha); Mxs(j,:,i) = r_curve*Q_force(j,:,i).*sin(alpha).*sin(phi); Mys(j,:,i) = -r_curve*Q_force(j,:,i).*sin(alpha).*cos(phi); % Stiffness Kxx = k*(A-A0).^n.*cos(phi).^2.*(n*A.*(delta_rs.^2)./(A-A0)+A.^2-delta_rs.^2)./(A.^3); kxx(j,i) = sum(Kxx(I)); Kxy = k*(A-A0).^n.*sin(phi).*cos(phi).*(n*A.*(delta_rs.^2)./(A-A0)+A.^2-delta_rs.^2)./(A.^3); kxy(j,i) = sum(Kxy(I)); Kxz = k*(A-A0).^n.*delta_rs.*delta_zs.*cos(phi).*(n*A./(A-A0)-1)./(A.^3); kxz(j,i) = sum(Kxz(I)); Kxrx = r_curve*k*(A-A0).^n.*delta_rs.*delta_zs.*sin(phi).*cos(phi).*(n*A./(A-A0)-1)./(A.^3); kxrx(j,i) = sum(Kxrx(I)); Kxry = r_curve*k*(A-A0).^n.*delta_rs.*delta_zs.*cos(phi).^2.*(1-n*A./(A-A0))./(A.^3); kxry(j,i) = sum(Kxry(I)); Kyy = k*(A-A0).^n.*sin(phi).^2.*(n*A.*(delta_rs.^2)./(A-A0)+A.^2-delta_rs.^2)./(A.^3); kyy(j,i) = sum(Kyy(I)); Kyz = k*(A-A0).^n.*delta_rs.*delta_zs.*sin(phi).*(n*A./(A-A0)-1)./(A.^3); kyz(j,i) = sum(Kyz(I)); Kyrx = r_curve*k*(A-A0).^n.*delta_rs.*delta_zs.*sin(phi).^2.*(n*A./(A-A0)-1)./(A.^3); kyrx(j,i) = sum(Kyrx(I)); Kyry = r_curve*k*(A-A0).^n.*delta_rs.*delta_zs.*sin(phi).*cos(phi).*(1-n*A./(A-A0))./(A.^3); kyry(j,i) = sum(Kyry(I)); Kzz = k*(A-A0).^n.*(n*A.*(delta_zs.^2)./(A-A0)+A.^2-delta_zs.^2)./(A.^3); kzz(j,i) = sum(Kzz(I)); Kzrx = r_curve*k*(A-A0).^n.*sin(phi).*(n*A.*(delta_zs.^2)./(A-A0)+A.^2-delta_zs.^2)./(A.^3); kzrx(j,i) = sum(Kzrx(I)); Kzry = r_curve*k*(A-A0).^n.*cos(phi).*(delta_zs.^2-n*A.*(delta_zs.^2)./(A-A0)-A.^2)./(A.^3); kzry(j,i) = sum(Kzry(I)); Krxrx = r_curve^2*k*(A-A0).^n.*sin(phi).^2.*(n*A.*(delta_zs.^2)./(A-A0)+A.^2-delta_zs.^2)./(A.^3); krxrx(j,i) = sum(Krxrx(I)); Krxry = r_curve^2*k*(A-A0).^n.*sin(phi).*cos(phi).*(delta_zs.^2-n*A.*(delta_zs.^2)./(A-A0)-A.^2)./(A.^3); krxry(j,i) = sum(Krxry(I)); Kryry = r_curve^2*k*(A-A0).^n.*cos(phi).^2.*(n*A.*(delta_zs.^2)./(A-A0)+A.^2-delta_zs.^2)./(A.^3); kryry(j,i) = sum(Kryry(I));