Используйте метод конечных элементов для решения двумерного уравнения диффузии (уравнение теплопроводности), но расчлените

Я пытаюсь использовать конечный элемент для решения двумерного уравнения диффузии:

numx = 101;   % number of grid points in x
numy = 101; 
numt = 1001;  % number of time steps to be iterated over 
dx = 1/(numx - 1);
dy = 1/(numy - 1);
dt = 1/(10*(numt - 1));


x = 0:dx:1;   %vector of x values, to be used for plotting
y = 0:dy:1;   %vector of y values, to be used for plotting

P = zeros(numx,numy,numt);   %initialize everything to zero

% Initial Condition
mu = 0.5;
sigma = 0.05;
for i=2:numx-1
    for j=2:numy-1
        P(i,j,1) = exp(-( (x(i)-mu)^2/(2*sigma^2) + (y(j)-mu)^2/(2*sigma^2))/(2*sigma^2)) / sqrt(2*pi*sigma^2);           
 end

 end

% Diffusion Equation
for k=1:numt
   for i=2:numx-1
       for j=2:numy-1
           P(i,j,k+1) = P(i,j,k) + .5*(dt/dx^2)*( P(i+1,j,k) - 2*P(i,j,k) + P(i-1,j,k) ) + .5*(dt/dy^2)*( P(i,j+1,k) - 2*P(i,j,k) + P(i,j-1,k) ); 
       end
   end
end

figure(1)
surf(P(:,:,30))

Начальное условие P распределено по Гауссу.

Однако при увеличении k (k=30 в моем коде) P взрывается. Однако P должно уменьшаться, так как это решение уравнения диффузии.

Я понятия не имею о проблеме в P(i,j,k+1), как решить эту проблему?

Спасибо


matlab differential-equations finite-element-analysis
person sleeve chen    schedule 25.10.2017    source источник
comment
Эта задача может иметь решение в закрытой форме. Если это не сработает, это может быть признаком того, что у вашего решения FEA нет шансов на успех. Убедитесь, что граничные и начальные условия имеют смысл.   -  person duffymo    schedule 14.11.2019

Ответы (1)


arrow_upward
0
arrow_downward

Значение P(1, :, :), P(101, :, :), P(:, 1, :) и P(:, 101 , :) в вашем коде всегда равно нулю. Это означает, что P всегда равен нулю на границе двумерного прямоугольника. Он не подчиняется закону диффузии.

Во-первых, вы должны задать начальное условие для границы.

% Initial Condition
mu = 0.5;
sigma = 0.05;
for i=1:numx
    for j=1:numy
        P(i,j,1) = exp(-( (x(i)-mu)^2/(2*sigma^2) + (y(j)-mu)^2/(2*sigma^2))/(2*sigma^2)) / sqrt(2*pi*sigma^2);           
     end
 end

А затем рассчитать значение Р на границе в процессе диффузии.

% Diffusion Equation
for k=1:numt
   for i=2:numx-1
       for j=2:numy-1
           P(i,j,k+1) = P(i,j,k) + .5*(dt/dx^2)*( P(i+1,j,k) - 2*P(i,j,k) + P(i-1,j,k) ) + .5*(dt/dy^2)*( P(i,j+1,k) - 2*P(i,j,k) + P(i,j-1,k) ); 
       end
   end

   % calculate the value of P at the boudary
   P(1, :, k + 1) = 2 * P(2, :, k + 1) - P(3, : , K + 1);
   P(numx, :, k + 1) = 2 * P(numx - 1, :, k + 1) - P(numx - 2, : , K + 1);
   P(:, 1, k + 1) = 2 * P(:, 2, k + 1) - P(:, 3 , K + 1);
   P(:, numy, k + 1) = 2 * P(:, numy - 1, k + 1) - P(:, numy - 2 , K + 1);
end
person ZHANG Luyao    schedule 30.10.2019