零通量边界下2D扩散模拟物质量不守恒的MATLAB技术问询

我使用MATLAB File Exchange上的2D扩散问题有限差分代码，设置所有边界条件为零通量，初始条件为u(20,25)=2，其余位置为0。运行约100个时间步后，预期每个时间步的总物质量（通过sum(sum(u))计算）保持为2，但结果显示总物质量随时间下降。请问是代码使用有误，还是对扩散概念的理解存在偏差？

% Simulating the 2-D Diffusion equation by the Finite Difference Method 
% Numerical scheme used is a first order upwind in time and a second order central difference in space (Implicit and Explicit)
%%
%Specifying parameters
nx=200;                           %Number of steps in space(x)
ny=200;                           %Number of steps in space(y)       
nt=100;                           %Number of time steps 
dt=0.01;                         %Width of each time step
dx=2/(nx-1);                     %Width of space step(x)
dy=2/(ny-1);                     %Width of space step(y)
x=0:dx:2;                        %Range of x(0,2) and specifying the grid points
y=0:dy:2;                        %Range of y(0,2) and specifying the grid points
u=zeros(nx,ny);                  %Preallocating u
un=zeros(nx,ny);                 %Preallocating un
vis=0.1;                         %Diffusion coefficient/viscocity
UW=0;                            %x=0 Dirichlet B.C 
UE=0;                            %x=L Dirichlet B.C 
US=0;                            %y=0 Dirichlet B.C 
UN=0;                            %y=L Dirichlet B.C 
UnW=0;                           %x=0 Neumann B.C (du/dn=UnW)
UnE=0;                           %x=L Neumann B.C (du/dn=UnE)
UnS=0;                           %y=0 Neumann B.C (du/dn=UnS)
UnN=0;                           %y=L Neumann B.C (du/dn=UnN)
%%
%Initial Conditions
u(nx/2,ny/2)=2;
%%
%B.C vector
bc=zeros(nx-2,ny-2);
bc(1,:)=UW/dx^2; bc(nx-2,:)=UE/dx^2;  %Dirichlet B.Cs
bc(:,1)=US/dy^2; bc(:,ny-2)=UN/dy^2;  %Dirichlet B.Cs
%bc(1,:)=-UnW/dx; bc(nx-2,:)=UnE/dx;  %Neumann B.Cs
%bc(:,1)=-UnS/dy; bc(:,nx-2)=UnN/dy;  %Neumann B.Cs
%B.Cs at the corners:
bc(1,1)=UW/dx^2+US/dy^2; bc(nx-2,1)=UE/dx^2+US/dy^2;
bc(1,ny-2)=UW/dx^2+UN/dy^2; bc(nx-2,ny-2)=UE/dx^2+UN/dy^2;
bc=vis*dt*bc;
%Calculating the coefficient matrix for the implicit scheme
Ex=sparse(2:nx-2,1:nx-3,1,nx-2,nx-2);
Ax=Ex+Ex'-2*speye(nx-2);        %Dirichlet B.Cs
%Ax(1,1)=-1; Ax(nx-2,nx-2)=-1;  %Neumann B.Cs
Ey=sparse(2:ny-2,1:ny-3,1,ny-2,ny-2);
Ay=Ey+Ey'-2*speye(ny-2);        %Dirichlet B.Cs
%Ay(1,1)=-1; Ay(ny-2,ny-2)=-1;  %Neumann B.Cs
A=kron(Ay/dy^2,speye(nx-2))+kron(speye(ny-2),Ax/dx^2);
D=speye((nx-2)*(ny-2))-vis*dt*A;
%%
%Calculating the field variable for each time step
i=2:nx-1;
j=2:ny-1;

totalField=[];
totalField(1)=sum(sum(u));
for it=0:nt
    un=u;
    h=surf(x,y,u','EdgeColor','none');       %plotting the field variable
    shading interp
    axis ([0 2 0 2 0 2])
    title({['2-D Diffusion with {\nu} = ',num2str(vis)];['time (\itt) = ',num2str(it*dt)]})
    xlabel('Spatial co-ordinate (x) \rightarrow')
    ylabel('{\leftarrow} Spatial co-ordinate (y)')
    zlabel('Transport property profile (u) \rightarrow')
    drawnow; 
    refreshdata(h)
    
    %Uncomment as necessary
    %Implicit method:
    U=un;U(1,:)=[];U(end,:)=[];U(:,1)=[];U(:,end)=[];
    U=reshape(U+bc,[],1);
    U=D\U;
    U=reshape(U,nx-2,ny-2);
    u(2:nx-1,2:ny-1)=U;
    %Boundary conditions
    %Dirichlet:
    u(1,:)=UW;
    u(nx,:)=UE;
    u(:,1)=US;
    u(:,ny)=UN;
 
    totalField=[totalField sum(sum(u))];
end
plot(totalField)

你的代码使用有误，核心问题是你想设置零通量（Neumann）边界条件，但实际代码中启用的是Dirichlet边界条件，并且没有正确配置Neumann对应的矩阵和边界向量。

具体问题点： #

1. 边界条件配置错误：

   • 你注释了Neumann边界条件的bc设置代码，当前启用的是Dirichlet边界（u=0），这会导致边界处的物质量被强制设为0，自然造成总物质量下降。
   • 同时，系数矩阵Ax和Ay也使用的是Dirichlet版本，没有切换到Neumann对应的矩阵修改（注释掉的Ax(1,1)=-1; Ax(nx-2,nx-2)=-1;部分）。

2. 初始条件设置偏差：
   你描述的初始条件是u(20,25)=2，但代码里写的是u(nx/2,ny/2)=2（nx=200时对应(100,100)），不过这不是总物质量下降的原因。

修正步骤： #

• 切换到Neumann边界条件配置：
  1. 注释掉Dirichlet的bc设置代码，取消注释Neumann的bc设置：

     % bc(1,:)=UW/dx^2; bc(nx-2,:)=UE/dx^2;  %Dirichlet B.Cs
     % bc(:,1)=US/dy^2; bc(:,ny-2)=UN/dy^2;  %Dirichlet B.Cs
     bc(1,:)=-UnW/dx; bc(nx-2,:)=UnE/dx;  %Neumann B.Cs
     bc(:,1)=-UnS/dy; bc(:,ny-2)=UnN/dy;  %Neumann B.Cs
     

  2. 注释掉Dirichlet的系数矩阵代码，取消注释Neumann的矩阵修改：

     % Ax=Ex+Ex'-2*speye(nx-2);        %Dirichlet B.Cs
     Ax=Ex+Ex'-2*speye(nx-2);
     Ax(1,1)=-1; Ax(nx-2,nx-2)=-1;  %Neumann B.Cs
     % Ay=Ey+Ey'-2*speye(ny-2);        %Dirichlet B.Cs
     Ay=Ey+Ey'-2*speye(ny-2);
     Ay(1,1)=-1; Ay(ny-2,ny-2)=-1;  %Neumann B.Cs
     

  3. 注释掉Dirichlet边界赋值代码，因为Neumann边界条件已经通过矩阵和bc向量实现，不需要再强制设置边界值：

     % %Boundary conditions
     % %Dirichlet:
     % u(1,:)=UW;
     % u(nx,:)=UE;
     % u(:,1)=US;
     % u(:,ny)=UN;
     

• 修正初始条件（可选，根据你的需求）：
  将u(nx/2,ny/2)=2;改为u(20,25)=2;。

修正后，零通量边界条件会被正确应用，总物质量将保持守恒，不会随时间下降。

内容的提问来源于stack exchange，提问作者Cnine