轴承油膜温度场程序

主要应用超松弛方法求解轴承底面油膜温度场，通过MATLAB进行数值模拟。

% Twodimensional heat conduction
% Finite Volume Method
% SOR
clear all;
x=[];y=[];T=[];Told=[];Su=[];Sp=[];ap=[];ae=[];aw=[];as=[];an=[];
great = 1.e20;
lambda = 10; % thermal conductivity
alfa = 10; % heat transfer coefficient
dt = great; % Time step. If great stedy state
density = 6000;% density
cp = 500;% heat capacity
Lx = 0.12; % length x-direction
Ly = 0.12; % length y -direction
Tfluid = 20; % Fluid temperature
Tinit = 50; % Initial guess and top- and bottom tempearature
%cv_x = input（‘Number of CVs in x-direction = ‘）
%cv_y = input（‘Number of CVs in y-direction = ‘）
cv_x=10;cv_y=10;
ni = cv_x+2; % grid points x-direction
nj = cv_y+2; % grid points y-direction
dx = Lx/cv_x;
dy = Ly/cv_y;
x（1） = 0;
x（2）=dx/2;
for i = 3:ni-1
   x（i）=x（i-1）+dx;
end;
x（ni）=Lx;
y（1） = 0;
y（2）=dy/2;
for j = 3:nj-1
   y（j）=y（j-1）+dy;
end
y（nj）=Ly;
% Initial values and coefficients
for i = 1:ni
  for j = 1:nj
    T（ij） = Tinit;  %Initial temperature
    Told（ij） = Tinit;
    T（i1） = 50;
    T（inj） = 50;
    Su（ij）=0;       %Initial indendendent source term
    Sp（ij）=0;       %Initial dependent source term
    ae（ij） = lambda*dy/dx;
    aw（ij） = lambda*dy/dx;
    an（ij） = lambda*dx/dy;
    as（ij） = lambda*dx/dy;
    dV = dx*dy;
    ap0 = density*cp*dV/dt;
    if i==2  % convective heat transfer boundary
       Su（ij） = Tfluid/（1/alfa+dx/（2*lambda））*dy/dV;
       Sp（ij） = -1/（1/alfa+dx/（2*lambda））*dy/dV;
       aw（ij） = 0;
    end;
    if i==ni-1 % insulated boundary
       ae（ij） = 0;
    end
    if j==2 % bottom boundary given temperature
       as（ij）=2*lambda*dx/dy;
    end   
    if j==nj-1 % top boundary given temperature
       an（ij）=2*lambda*dx/dy;
    end   
    ap（ij） = ae（ij）+aw（ij）+an（ij）+as（ij）-Sp（ij）*dV+ap0;
  end;
end;
%%%%%%%%%%%
maxres = 1.0e-6;
maxit = 100;
time=0;
maxtime=100;
s=（cos（pi/cv_x）+（dx/dy）^2*cos（pi/cv_y））/（1+（dx/dy）^2）;
omega =2/（1+sqrt（1-s^2））;omega=1;
while （time < （maxtime+dt/2））
 Told=T;  
 sumres = 1;
 counter = 0;
 while （sumres>maxres&counter  sumres = 0;
  for i = 2:ni-1
   for j = 2:nj-1
    T（ij）=omega*（ae（ij）*T（i+1j）+aw（ij）*T（i-1j）+an（ij）*T（ij+1）...
       +as（ij）*T（ij-1）+Su（ij）*dV+ap0*Told（ij））/ap（ij）+（1-omega）*T（ij）;
    res =  abs（ap（ij）*T（ij）-（ae（ij）*T（i+1j）+aw（ij）*T（i-1j）+...
      an（ij）*T（ij+1）+as（ij）*T（ij-1）+Su（ij）*dV+ap0*Told（ij）））;
    sumres=sumres+res;
   end;
  end;
  for i = 2:ni-1
   for j = 2:nj-1
    res =  abs（ap（ij）*T（ij）-（ae（ij）*T（i+1j）+aw（ij）*T（i-1j）+...
      an（ij）*T（ij+1）+as（ij）*T（ij-1）+Su（ij）*dV+ap0*Told（ij）））;
    sumres=sumres+res;
   end;
  end

  sumerr=sumres
  counter = counter + 1
 end;
 time = time +dt;
end; 
% Calculate boundary values
 for j = 2:nj-1
   T（1j）=（alfa*Tfluid+lambda/（dx/2）*T（2j））/（alfa+lambda/（dx/2））;
   T（nij） = T（ni-1j）;
end;   
%
pcolor（xyT‘）;shading interp;xlabel（‘x‘）;ylabel（‘y‘）;title（‘Tempera

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3036  2003-08-07 14:21  temperatureofoilthrustbearing\Example_3_SOR.m

     目录          0  2008-11-07 12:16  temperatureofoilthrustbearing

----------- ---------  ---------- -----  ----

                 3254                    3