基于matlab有限元法实现动压润滑推力轴承的简单分析附代码-CSDN博客

本文链接：https://blog.csdn.net/Matlab_dashi/article/details/135900804

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，

代码获取、论文复现及科研仿真合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab完整代码及仿真定制内容点击👇

智能优化算法神经网络预测雷达通信无线传感器电力系统

信号处理图像处理路径规划元胞自动机无人机

🔥 内容介绍

摘要

本文提出了一种基于有限元法（FEM）的动压润滑推力轴承简单分析方法。该方法将推力轴承离散为有限个单元，并应用适当的边界条件和载荷条件，然后求解单元内的控制方程组，得到每个单元的压力分布和位移场。最后，通过积分得到轴承的总载荷和摩擦力矩。该方法可以用于分析不同几何形状、材料和工况条件下的推力轴承的性能。

引言

推力轴承是一种广泛应用于旋转机械中的轴承，其主要作用是承受轴向载荷。动压润滑推力轴承是利用流体动压效应来实现润滑的，具有摩擦力矩小、承载能力大、寿命长等优点。因此，动压润滑推力轴承在航空发动机、燃气轮机、压缩机等高速旋转机械中得到了广泛的应用。

有限元法简介

有限元法（FEM）是一种广泛应用于工程分析中的数值方法。其基本思想是将连续介质离散为有限个单元，并应用适当的边界条件和载荷条件，然后求解单元内的控制方程组，得到每个单元的压力分布和位移场。最后，通过积分得到整个结构的总载荷和位移。

动压润滑推力轴承的有限元分析

在动压润滑推力轴承的有限元分析中，首先需要将推力轴承离散为有限个单元。单元的形状和大小可以根据推力轴承的几何形状和载荷分布来确定。常用的单元形状有三角形、四边形和六边形。

然后，需要为每个单元应用适当的边界条件和载荷条件。边界条件可以是位移边界条件、力边界条件或混合边界条件。载荷条件可以是集中载荷、分布载荷或混合载荷。

接下来，需要求解单元内的控制方程组。控制方程组通常包括连续性方程、动量方程和能量方程。求解控制方程组的方法有很多种，常用的方法有直接法和迭代法。

最后，通过积分得到轴承的总载荷和摩擦力矩。总载荷是轴承承受的所有载荷的总和，摩擦力矩是轴承旋转时所产生的摩擦阻力。

📣 部分代码

function [p,pz,A,DIA_IN,DIA_OUT,idx_boundary_in,idx_boundary_out,fx,fz,Kzz] = ReyThrustStiffFunc(...PAD_DIM,TH_DIM,RA_DIM,ns,es,Dt,Dr,HP,AS,VISCO,VIS_EN,TURB_SWITCH,RHO,PB,ALPHA,DISP_PREFIX,ANG_OFF)% 调试开关IS_DEBUG = 1;% 计算常量ELE_NODES_NUM = 4;        % 每个单元的节点个数，四节点STIFFNESS_BASE = 1e9;     % 用于将求得的刚度值以10^9形式显示% DAMPING_BASE   = 1e6;     % 用于将求得的阻尼值以10^6形式显示% 求解一次稳态压力场STABLE_FIELD_DISP_PREFIX = [DISP_PREFIX,'Stable pressure field: '];[p,A,DIA_IN,DIA_OUT,idx_boundary_in,idx_boundary_out,fx,fz] = ReyThrustFunc(...PAD_DIM,TH_DIM,RA_DIM,ns,es,Dt,Dr,HP,AS,VISCO,VIS_EN,TURB_SWITCH,RHO,PB,ALPHA,0,STABLE_FIELD_DISP_PREFIX,ANG_OFF);% ========================= 组装方程右端 ==============================% 转换有限元单元积分参数Dtn = 1/Dt;Drn = 1/Dr;det_J = Dt*Dr;% 用于表示方向的变量a   = -1;atd = 1;tRe = RHO*AS/(VISCO*VIS_EN); % 湍流计算时局部雷诺数计算临时变量% 方程右端临时向量，列向量rhs_idx = zeros(1,4); % 方程右端索引% 周向（第一坐标）和径向（第二坐标）的方程RHS临时变量tRHSAr = -1*a/(VISCO*VIS_EN) * Drn^2 * det_J * [1/3,-1/3,1/3,-1/3]';tRHSAt = -1*a/(VISCO*VIS_EN) * Dtn^2 * det_J * [1/3,-1/3,1/3,-1/3]';tRHSB  = -1*a/(VISCO*VIS_EN) * Dtn^2 * det_J * [-1,1,1,-1]'; % t方向tRHSC  = -1*a/(VISCO*VIS_EN) * Drn^2 * det_J * [-1,-1,1,1]'; % r方向% 阻尼计算的临时方程右端% tRHSzt = atd * det_J * [1,1,1,1]';% RHSRHSz  = zeros((TH_DIM+1)*(RA_DIM+1),1);% RHSzt = RHSz;% 止推瓦块的起始ths = PAD_DIM(1,1);thend = PAD_DIM(1,2);% 止推轴承湍流系数计算用的参数MS  = -0.25;NS  = 0.066;tGx = 2^(1+MS) / (NS * (2 + MS));tGz = 2^(1+MS) / NS;if(IS_DEBUG)    check_reynolds = zeros(1,TH_DIM*RA_DIM); % 记录局部雷诺数用于调试    check_rac      = zeros(1,TH_DIM*RA_DIM); % 记录单元中心位置用于调试    check_h        = zeros(1,TH_DIM*RA_DIM); % 记录液膜厚度用于调试end% 显示信息disp([DISP_PREFIX,'Assemble RHSs of perturbation equations...']);for I = 1:1:TH_DIM*RA_DIM    % 组装每一个单元        % 取出四个节点    for J = 1:1:ELE_NODES_NUM        tn(J,:) = [es(I,J),ns(es(I,J),1),ns(es(I,J),2)];    end % J        % 计算单元的中心坐标, theta center和axis center    thc = (tn(1,2) + tn(2,2))/2 + ANG_OFF;    rac = (tn(1,3) + tn(4,3))/2;        % 计算液膜厚度    h = HP + ALPHA*rac*sin(thend - thc);        % 计算圆周方向的系数Gx和Gz    % 径向轴承的描述%     if(TURB_SWITCH == 0)%         Gx = 1 / (4 * rac);%         Gz = rac / 4;%     elseif(TURB_SWITCH == 1)%         lRe = tRe * h * rac;%         if(lRe > 2000)%             Gx = 1 / ((12+0.0136*lRe^0.9) * rac / 3);%             Gz = rac / ((12+0.0043*lRe^0.96) / 3);%         else%             Gx = 1 / (4 * rac);%             Gz = rac / 4; %         end % lRe > 2000%     end    % ===== 2012.3.31 =========    % 止推轴承的描述    if(TURB_SWITCH == 0)        Gx = 1 / (4 * rac);        Gz = rac / 4;    elseif(TURB_SWITCH == 1)        lRe = tRe * h * rac;        if(IS_DEBUG)            check_reynolds(1,I) = lRe;            check_rac(1,I) = rac;            check_h(1,I) = h;        end        gRe = lRe^(1+MS);        if(lRe > 2000)%             Gx = 1 / ((12+0.0136*lRe^0.9) * rac / 3);            Gx = tGx / gRe / rac * 3;%             Gz = rac / ((12+0.0043*lRe^0.96) / 3);            Gz = tGz / gRe * rac * 3;        else            Gx = 1 / (4 * rac);            Gz = rac / 4;         end % lRe > 2000    end        % 临时节点压力    p0t = p(tn(:,1),1);        p0A = sum(0.25 * [1,-1,1,-1]' .* p0t);    p0B = sum(0.25 * [-1,1,1,-1]' .* p0t);    p0C = sum(0.25 * [-1,-1,1,1]' .* p0t);        ttRHS = (tRHSAt .* p0A .* h^2 + tRHSB .* p0B .* h^2) * Gx;    trRHS = (tRHSAr .* p0A .* h^2 + tRHSC .* p0C .* h^2) * Gz;      % 计算方程右端    rhs_idx = [tn(1,1),tn(2,1),tn(3,1),tn(4,1)];    RHSz(rhs_idx,1) = RHSz(rhs_idx,1) + trRHS + ttRHS;%     RHSzt(rhs_idx,1) = RHSzt(rhs_idx,1) + tRHSzt * rac;end% 处理方程右端，特别注意边界！！！还未修改RHSz( idx_boundary_in) = PB*DIA_IN;RHSz(idx_boundary_out) = PB*DIA_OUT;% RHSzt( idx_boundary_in) = PB*DIA_IN;% RHSzt(idx_boundary_out) = PB*DIA_OUT;% 显示信息disp([DISP_PREFIX,'RHSs assembled.']);% =========================== Solve ==================================% % z方向的速度摄动压力场% pzt = A\RHSzt;% disp('pzt solved.');% z方向的位移摄动压力场pz  = A\RHSz;disp([DISP_PREFIX,'pz solved.']);% 查看结果res_vec = A*pz - RHSz;res_norm = norm(res_vec);disp([DISP_PREFIX,'The norm of residual vector of pz is ',num2str(res_norm),'.']);% res_vec = A*pzt - RHSzt;% res_norm = norm(res_vec);% % disp(['The norm of residual vector of pzt is ',num2str(res_norm),'.']);if(IS_DEBUG)    % 检查局部雷诺数的值    disp([DISP_PREFIX,'The maximun Reynolds number is ',...        num2str(max(check_reynolds)),'.']);        disp([DISP_PREFIX,'The maximun r number is ',...        num2str(max(check_rac)),'.']);        disp([DISP_PREFIX,'The maximun h number is ',...        num2str(max(check_h)),'.']);    % 清理临时变量    clear check_reynolds    clear check_rac    clear check_hend % IS_DEBUG% ============================= 显示结果 ==============================% 重构结果向量形成一个矩阵矩阵的为（TH_DIM+1）行，（RA_DIM+1）列pzr  = reshape(pz,TH_DIM+1,RA_DIM+1);  % 位移摄动压力场重构% pztr = reshape(pzt,TH_DIM+1,RA_DIM+1); % 速度摄动压力场重构% th_idx_re = ((1:1:(TH_DIM+1))-1) .* Dt ./ (2*pi) .* 360; % re for reconstruction% ra_idx_re = ((1:1:(RA_DIM+1))-1) .* Dr + PAD_DIM(1,3);     % re for reconstruction% % disp([DISP_PREFIX,'Display the results in 3D.']);% % figure% % surf(ra_idx_re,th_idx_re,pzr,'LineStyle','none')% figure% % surf(ra_idx_re,th_idx_re,pztr,'LineStyle','none')% ========================= 对摄动压力场求积分 ===========================% 梯形积分Kzz = 0; % 与载荷平行的刚度% Czz = 0; % 与载荷平行的阻尼Kzz_temp = 0; % 临时变量% Czz_temp = 0; % 临时变量% fth_begin = -1*Dt/2; % 临时变量，force theta% fth = 0; % 临时变量，角度位置，单位为radfa = PAD_DIM(1,3);for I = 1:1:RA_DIM%     fth = fth_begin;    for J = 1:1:TH_DIM%         fth = fth + Dt;        %         if(J ~= TH_DIM)        if(1)            % 未到回环边界            Kzz_temp = pzr(  J,I)*fa + pzr(  J,I+1)*(fa+Dr) + ...                       pzr(J+1,I)*fa + pzr(J+1,I+1)*(fa+Dr);%             Czz_temp = pztr(  J,I) + pztr(  J,I+1) + ...%                        pztr(J+1,I) + pztr(J+1,I+1);        else            % 到达回环边界            Kzz_temp = pzr(J,I) + pzr(J,I+1) + ...                       pzr(1,I) + pzr(1,I+1);%             Czz_temp = pztr(J,I) + pztr(J,I+1) + ...%                        pztr(1,I) + pztr(1,I+1);        end               Kzz = Kzz + (-1)*Kzz_temp;%        Czz = Czz + (-1)*Czz_temp;    end % J    fa = fa + Dr;end % I% 乘以积分系数Kzz = Kzz * (Dt * Dr)/4;% Czz = Czz * (Dt * Dr)/4;% 显示刚度信息disp([DISP_PREFIX,'Kzz = ',num2str(Kzz/STIFFNESS_BASE),'x10^9 N/m']);% disp(['Czz = ',num2str(Czz/DAMPING_BASE),'x10^6 Ns/m']);if(IS_DEBUG)    disp([DISP_PREFIX,'Debug mode is on.']);end