уравнение Лапласа для моделирования диффузии тепла через L-образный канал - PullRequest
Новая
       62

уравнение Лапласа для моделирования диффузии тепла через L-образный канал

0 голосов
/ 03 мая 2020

Привет, я пытаюсь решить уравнение Лапласа для следующей формы. [! [Введите описание изображения здесь] [1]] [1]

[1]: https://i.stack.imgur.com/HEyG3.png. Я испытываю трудности с реализацией матрицы A и b в для l oop в следующем коде Matlab. Матрица b представляет собой граничное условие, а матрица A представляет собой коэффициенты уравнения с использованием метода центральной разности. 

    % Problem 2
clear;
% Set h = dx = dy
h = 0.025; %Try different values of h, until observing convergence.

% Geometry & boundary conditions
xmin = 0.0; x1 = 1.0; xmax = 2.0;
ymin = 0.0; y1 = 1.0; ymax = 2.0;
% Boundary conditions (In this problem, all boundaries are isothermal)
Tin = 100; %top
Tout = 0; %right
% Mesh (h is given as input)
i1 = (x1-xmin)/h+1; % i-index of nodes with x = x1;
j1 = (y1-ymin)/h+1; % j-index of nodes with y = y1;
imax = (xmax-xmin)/h; % i-index of the (interior) node with largest x.
jmax = (ymax-ymin)/h; % j-index of the (interior) node with largest y.
N = i1*j1 + (imax-i1)*j1 + i1*(jmax-j1); % total number of nodes in domain interior & on adiabatic walls
% Allocate memory for A and b (sparse matrix & vector).
A = sparse(N,N); b = sparse(N,1);
% Loop through all the interior nodes, and fill A & b.
for k=1:N
    A(k,k) = 4; % the diagonal element, corresponding to Tij, is always 4
    %In the following we look at the four neighbours of (i,j) that are
    %involved in the 5-point formula
    % Get (i,j) from k
    if (k<=imax*j1) %y<=1.0
        i = mod(k-1,imax) + 1;
        j = (k-i)/imax + 1;
    else %1.0<y<=2.0
        i = mod(k-imax*j1-1, i1) + 1;
        j = j1 + (k-imax*j1-i)/i1 + 1;
    end
A=?
B=?
end

% Solve Ax = b for nodal temperature values
T = full(A\b);

% Output Temperature at (1.0, 0.5)
i = (0.5*(xmax-xmin))/h+1;
j = (0.25*(ymax-ymin))/h+1;
Tsensor = T(imax*(j-1)+i);
fprintf('h = %e, T(%d,%d) = %e.\n', h, i, j, Tsensor);

% Visualization
figure
[X, Y] = meshgrid(xmin:h:xmax, ymin:h:ymax);
Tvis = zeros(size(X)); %just for visualization
Tvis(:,:) = NaN;
for j=1:j1
  for i=1:imax
      Tvis(i,j) = T(imax*(j-1)+i);
  end
  Tvis(imax+1,j) = Tout;
end
for j=j1+1:jmax
  for i=1:i1
      Tvis(i,j) = T(imax*i1+(j-j1-1)*i1+i);
  end
end
Tvis(1:i1,jmax+1) = Tin;

h = surf(X,Y,Tvis');
colormap('jet');
view(0,90);
shading interp
set(h,'edgecolor','k');
cb = colorbar; cb.Label.String = 'Temperature (C)';
xlabel('X (m)')
ylabel('Y (m)')
zlabel('Temperature (C)')
daspect([1 1 1])
...