我是 CFD 和单相/多相流和运输的新手。作为我学习探索的一部分，我正在尝试对二维域中心的固体物体周围的溶质传输进行建模。控制方程是对流-扩散方程，给出如下：
其中 C 是浓度，v 是流速，D 是扩散系数。在恒定流速下，上述简化为：
到目前为止，这是我的代码：
clc; clear; close all;
%%%%%%%%%% Specify inputs
CircleDiam = 40; %in pixels
Co = 0; %initial conc. in the domain [mol/L]
Cin = 0.01; %conc. of injected fluid [mol/L]
Lx = 0.1/100; %Length of domain [m]
nx = 200; %spatial gridpoints in x
dx = Lx/(nx-1); %Length step size [m]
Ly = 0.05/100; %Width of domain [m]
ny = 100; %spatial gridpoints in y
dy = Ly/(ny-1); %Length step size [m]
T = 5/(3600*24); %Simulation time [days]
nt = 8000; %shifts
dt = T/nt; %Time step [days]
% Flow
u = 103.68; %Velocity in x direction [m/day]
v = 0; %Velocity in y direction [m/day]
De = 8.64e-04; %Dispersion coeff. [m2/day]
betaX = u*dt/dx;
betaY = v*dt/dy;
gammaX = De*dt/(dx^2);
gammaY = De*dt/(dy^2);
%%%%%%%%%% Create image with solid object
radius = CircleDiam/2;
% obtain full output grids from grid vectors
[Colgrids, Rowgrids] = meshgrid(1:nx, 1:ny);
% create a logical mask for the circle by specifying the center and diameter of the circle.
centerX = 0.5 + (nx/2);
centerY = 0.5 + (ny/2);
% obtain image from: ( (y-y0)^2 + (x-x0)^2 ) <= r^2, where (y0,x0) is the centre point of circle
SolidImg = (Rowgrids - centerY).^2 + (Colgrids - centerX).^2 <= radius.^2;
% change from logical to numeric labels. Also, transpose matrix to conform with the conc. matrix
P = double(SolidImg');
%figure, imshow(~SolidImg, [], 'InitialMagnification','fit'); box on;
% Gridblocks
x = 0:dx:Lx;
y = 0:dy:Ly;
t = 0:dt:T;
[X,Y] = meshgrid(x, y);
%specify initial conditions
C = zeros(nx, ny, nt+1);
C(:,:,1) = Co; %Initial condition
%iterate finite difference equations
for k = 1:nt
for j = 2:ny-1
for i = 2:nx-1
if P(i,j)==1 %Solid pixels
%C(i,j,k+1) = C(i,j,k);
C(i,j,k+1) = 0;
else
C(i,j,k+1) = C(i-1,j,k)*(betaX/2+gammaX) + C(i+1,j,k)*(gammaX-betaX/2)...
+ C(i,j-1,k)*(betaY/2+gammaY) + C(i,j+1,k)*(gammaY-betaY/2)...
+ C(i,j,k)*(1-2*gammaX-2*gammaY);
end
end
end
% Insert boundary conditions
C(:,1,k+1) = C(:,2,k+1); % bottom bc
C(:,end,k+1) = C(:,end-1,k+1); % top bc
C(1,:,k+1) = Cin; % left bc
C(end,:,k+1) = C(end-1,:,k+1); % right bc
% C_all = squeeze(C(:,:,k));
% pcolor(C_all'); shading interp; colorbar; colormap(jet); hc = colorbar;
% ylabel(hc,'Concentration [mol/L]');
% caxis([0 max(max(max(C)))])
% title( sprintf('Time = %f seconds', k) )
% pause(0.001);
%
end
我的结果（见第一张图片）看起来不像我期望的那样。理想情况下，我希望看到类似上面（第二个）图像的东西，所以我想知道我可能在哪里弄错了。我假设在固体边界处速度和扩散系数都为零。在实体边界处是否需要实施额外的边界条件？
如果有人可以在这里帮助我或为我指明正确的方向，我将不胜感激。
PS：
- 两个图像中的轴相同。不同之处在于我的以像素为单位，而另一个以厘米为单位。
谢谢你的期待。
编辑：
- 以前，我有 x 方向的速度 = y 方向的速度。这是错误的，因为流量应该只来自一个边缘。我现在已将其更正为
velocity in x = +103.68 m/day, 和velocity in y = 0。 - 以前，底部边缘的边界条件是
C(:,2,k+1) = C(:,1,k+1);。它应该是C(:,1,k+1) = C(:,2,k+1);。 - 上述更改意味着我现在有一个比我以前的版本更接近问题的解决方案（见第一张图片）。