使用有限差分法建模三角形误差

我想根据其尺寸在正方形单元格中建模一个三角形。但是，我对代码的某些部分有一些问题，我认为这可能是我编程背后的逻辑，或者是代码本身的定义。使用有限差分法建模三角形误差

% Simulation for triangle 

clear all; 
clc 
clear window 

%% INITIAL PARAMETERS 

% Variables 
ni=21;    % number of nodes in the cell - x axis 
nj=21;    % number of nodes in the cell - y axis 
delta_sq=0.3;   % microasperity area density (MAXIMUM 0.433) 
u=3.5;     % slider speed [m/s] 
v=0.8;     % lubricant viscosity [Pas] 
F=98.1;     % applied load [N/m^2] 
h_init=20e-06;   % initial guess of film thickness [m] 
e_converge=0.0001;  % convergence limit 
i_max=500    % maximum number of iterations 

% Constants 
r_1=600e-06;   % length of half unit cell [m] 
P_out=0;    % outer pressure boundary [N/m^2] 
P_init=0;    % inner pressure boundary [N/m^2] 
P_cav=0;    % cavitation pressure [N/m^2] 


%% MESH 
R_X=2*r_1;    % length of unit cell - x axis 
R_Y=2*r_1;    % length of unit cell - y axis 
dX=R_X/(ni-1);   % number of mesh of unit cell - x axis 
dY=R_Y/(nj-1);   % number of mesh of unit cell - y axis 
X(1)=-R_X/2;    % film thickness at first node - x axis 
Y(1)=-R_Y/2;    % film thickness at first node - y axis 

for ii=2:(2*ni-1), 
    X(ii)=X(ii-1)+dX/2; % location of intermediate nodes - x axis 
end 

for jj=2:(2*nj-1), 
    Y(jj)=Y(jj-1)+dY/2; % location of intermediate nodes - y axis 
end 


%% GEOMETRICAL BOUNDARIES 

L_triangle=sqrt((delta_sq*16*(r_1)^2)/(sqrt(3)));% side length of triangle 

h_step = 5e-06; % step microasperity depth 

nHi=max(size(X)); 
nHj=max(size(Y)); 

% flat step condition 
for ii=1:nHi, 
    for jj=1:nHj, 
     h(ii,jj)=h_init; 
     if X(ii)>(-L_triangle/3) & X(ii)<0; 
      if Y(jj)>0; 
       if Y(jj)>0 & Y(jj)<((L_triangle+3*X(ii))/(sqrt(3))); 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
      if Y(jj)<0; 
       if Y(jj)>(-L_triangle/(2*sqrt(3))) & Y(jj)<0; 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
     end 
     if X(ii)>0 & X(ii)<(L_triangle/3); 
      if Y(jj)>0; 
       if Y(jj)>0 & Y(jj)<((L_triangle-3*X(ii))/sqrt(3)); 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
      if Y(jj)<0; 
       if Y(jj)>(-L_triangle/(2*sqrt(3))) & Y(jj)<0; 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
     end 
     if X(ii)>0 & X(ii)<(L_triangle/2) ; 
      if Y(jj)<0; 
       if Y(jj)>(((-3*X(ii))+L_triangle)/sqrt(3)) & Y(jj)<0; 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
     end 
     if X(ii)>(-L_triangle/2) & X(ii)<0; 
      if Y(jj)<0; 
       if Y(jj)>(-sqrt(3)*(-3*X(ii)-L_triangle)/3) & Y(jj)<0; 
        h(ii,jj)=h_init+h_step; 
       else 
        h(ii,jj)=h_init; 
       end 
      end 
     end 
    end 
end 






%% Pressure 

k=1;    
for ii=1:ni, 
    for jj=1:nj, 
     P(ii,jj,k)=P_init; 
end 
end 

%Define Coefficients 
i=0; 
for ii=1:2:nHi, 
    i=i+1; 
    j=1; 
    if ii==1, 
     for jj=3:2:nHj-2, 
     j=j+1; 
      A(i,j)=((6*v*u/dX)*(h(ii+1,jj)-h(nHi-1,jj)))/((h(nHi-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      B(i,j)=((h(ii+1,jj)^3)/dX^2)/((h(nHi-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      C(i,j)=((h(nHi-1,jj)^3)/dX^2)/((h(nHi-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      D(i,j)=((h(ii,jj+1)^3)/dY^2)/((h(nHi-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      E(i,j)=((h(ii,jj-1)^3)/dY^2)/((h(nHi-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
     end 
     elseif ii==nHi, 
     for jj=3:2:nHj-2, 
     j=j+1; 
      A(i,j)=((6*v*u/dX)*(h(2,jj)-h(ii-1,jj)))/((h(ii-1,jj)^3+h(2,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      B(i,j)=((h(2,jj)^3)/dX^2)/((h(ii-1,jj)^3+h(2,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      C(i,j)=((h(ii-1,jj)^3)/dX^2)/((h(ii-1,jj)^3+h(2,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      D(i,j)=((h(ii,jj+1)^3)/dY^2)/((h(ii-1,jj)^3+h(2,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      E(i,j)=((h(ii,jj-1)^3)/dY^2)/((h(ii-1,jj)^3+h(2,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
     end 
     else 
    for jj=3:2:nHj-2, 
     j=j+1; 
      A(i,j)=((6*v*u/dX)*(h(ii+1,jj)-h(ii-1,jj)))/((h(ii-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      B(i,j)=((h(ii+1,jj)^3)/dX^2)/((h(ii-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      C(i,j)=((h(ii-1,jj)^3)/dX^2)/((h(ii-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      D(i,j)=((h(ii,jj+1)^3)/dY^2)/((h(ii-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
      E(i,j)=((h(ii,jj-1)^3)/dY^2)/((h(ii-1,jj)^3+h(ii+1,jj)^3)/dX^2+(h(ii,jj-1)^3+h(ii,jj+1)^3)/dY^2); 
    end 
    end 
end 

%Solution Kernal 
err=1.0; 
while err>e_converge, 
%Set the boundary conditions 
    for ii=1:ni, 
    P(ii,1,k)=P_init; 
    P(ii,nj,k)=P_out; 
    end 
    for i=1:ni, 
    for j=2:nj-1, 
     if i==1, 
      P(i,j,k+1)=(A(i,j)+B(i,j)*P(i+1,j,k)+C(i,j)*P(ni,j,k)+D(i,j)*P(i,j+1,k)+E(i,j)*P(i,j-1,k)); 
      elseif i==ni, 
      P(i,j,k+1)=(A(i,j)+B(i,j)*P(1,j,k)+C(i,j)*P(i-1,j,k)+D(i,j)*P(i,j+1,k)+E(i,j)*P(i,j-1,k)); 
     else 
      P(i,j,k+1)=(-A(i,j)+B(i,j)*P(i+1,j,k)+C(i,j)*P(i-1,j,k+1)+D(i,j)*P(i,j+1,k)+E(i,j)*P(i,j-1,k+1)); 
     end 
      if P(i,j,k+1)<P_cav 
       P(i,j,k+1)=P_cav; 
      end 
    end 
    end 


% error check point 
    PP=max(max(P(:,:,k+1))); 
    LIM=0; 
    for i=1:ni, 
    for j=2:nj-1, 
     conv=(P(i,j,k+1)-P(i,j,k))/PP; 
     LIM=LIM+conv^2; 
    end 
    end 
    err=1/((ni)*(nj-2))*sqrt(LIM); 
    k=k+1; 
    if k>i_max, 
     k 
    break 
    end 
end 

% Final iteration 
for ii=1:ni, 
    P(ii,1,k)=P_init; 
    P(ii,nj,k)=P_out; 
end 
for ii=1:ni, 
    for jj=1:nj, 
     P_solution(ii,jj)=P(ii,jj,k); 
    end 
end 

i=0; 
for ii=1:2:nHi, 
    i=i+1; 
    x_i(i)=X(ii); 
    y_j(i)=Y(ii); 
end 

Maximum_Pressure=max(max(P_solution)) 
Average_Pressure=mean(mean(P_solution)) 
Minimum_Pressure=min(min(P_solution)) 

figure 
surf(X,Y,h) 

figure 
surf(y_j,x_i,P_solution)

我附上了代码。问题在于三角形的后半部分。我还附上了我分析的基础。问题在于代码的-s/2到s/2之间。我会赞赏任何关于我的代码或定义的错误逻辑的暗示。单元格是2r1 x 2r1尺寸的正方形尺寸

当您运行代码时，三角形的上半部分看起来不错，但下半部分没有。我附上了我的代码背后的分析图像。

http://postimg.org/image/kvn2xqlcx/

来源

2014-09-01 blueant1993

这是很多代码。你需要做的是仔细检查它（逐步或设置断点），并在此时开始偏离你的期望。 MATLAB有很多这样的内置工具：http://www.mathworks.co.uk/help/matlab/debugging-code.html – nkjt 2014-09-01 09:57:52

一般情况下，这种问题是气馁上堆栈溢出，因为这个问题方面的一个非常具体的代码调试，因此将用处不大将来其他人。

对此代码一些提示：

你不应该忽略lint警告（橙色强调你在MATLAB编辑器中看到） - 他们是有原因的。 There's a known joke on this topic...
如果您的代码不是clear all，而是使用clear variables--原因是clear all也会清除断点（您通常希望保留它）。
我建议将close all force添加到初始化，以删除任何以前打开的数字。
你应该仔细看看vectorization - 这可以（通常）帮助你的代码运行得更快，并且更具可读性。

话虽如此 - 这应该解决您的问题（以下的取代你行89-106）：

if X(ii)>0 && X(ii)<(L_triangle/2); 
     if Y(jj)<0; 
      if Y(jj)>(((-3*X(ii))+L_triangle)/sqrt(3)); 
       h(ii,jj)=h_init; 
      else 
       h(ii,jj)=h_init+h_step; 
      end 
     end 
    end 
    if X(ii)>(-L_triangle/2) && X(ii)<0; 
     if Y(jj)<0; 
      if Y(jj)>(-sqrt(3)*(-3*X(ii)-L_triangle)/3); 
       h(ii,jj)=h_init; 
      else 
       h(ii,jj)=h_init+h_step; 
      end 
     end 
    end 
    if X(jj)<(-L_triangle/(2*sqrt(3))) 
     h(ii,jj)=h_init; 
    end

你在错误的地点是在加入+h_step，也是区域“下的治疗“三角失踪了。

来源

2014-09-01 10:07:01

使用有限差分法建模三角形误差

回答

相关问题