Моделирование треугольника с использованием метода конечных разностей

Question

Я пытаюсь смоделировать треугольник в квадратной ячейке на основе его измерения. Тем не менее, у меня есть некоторые проблемы с частями кода, я думаю, что это может быть логикой моего программирования или над определением с самим кодом.

% 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

при запуске кода верхняя половина треугольника выглядит нормально, но нижняя половина isnt. Я прикрепил образ анализа моего кода.

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

Answer 1

Как правило, этот вид вопросов обескуражен на переполнение стека, так как вопрос касается отладки весьма специфического кода и, следовательно, будет мало пользы для других людей в будущем.

Некоторые советы относительно этого кода:

1. Вы не должны игнорировать предупреждения пуха (оранжевый подчеркивает вы видите в редакторе MATLAB) - они там по причине. There's a known joke on this topic...
2. Ваш код будет намного проще отлаживать, если вместо clear all вы используете clear variables - причина в том, что clear all также очищает контрольные точки (которые вы обычно хотите сохранить).
3. Я бы предложил добавить к инициализации close all force, чтобы удалить все ранее открытые фигуры.
4. Вы должны действительно взглянуть на 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 в неправильном месте, а также лечение области «под «Треугольник пропал.