用有限元法模拟地下水，matlab代码

以下是一个简单的用有限元法模拟地下水的matlab代码：

% 定义模拟区域
L = 100; % 模拟区域长度
H = 50; % 模拟区域宽度
nx = 10; % x方向网格数
ny = 5; % y方向网格数
dx = L/nx; % x方向网格大小
dy = H/ny; % y方向网格大小

% 定义材料参数
K = 10; % 渗透率
S = 0.1; % 孔隙度
phi = S; % 孔隙率
rho = 1000; % 密度
g = 9.81; % 重力加速度

% 定义时间步长和模拟时间
dt = 0.01; % 时间步长
T = 100; % 模拟总时间
nt = T/dt; % 时间步数

% 初始化水头和流量矩阵
h = zeros(nx*ny,1);
q = zeros(nx*ny,1);

% 构造有限元刚度矩阵和质量矩阵
A = zeros(nx*ny,nx*ny);
M = zeros(nx*ny,nx*ny);

for i = 1:nx
    for j = 1:ny
        n = (i-1)*ny+j;
        if i > 1
            A(n,n-ny) = A(n,n-ny) - K/dx;
            M(n,n-ny) = M(n,n-ny) + phi*rho*dx*dy/dt;
        end
        if i < nx
            A(n,n+ny) = A(n,n+ny) - K/dx;
            M(n,n+ny) = M(n,n+ny) + phi*rho*dx*dy/dt;
        end
        if j > 1
            A(n,n-1) = A(n,n-1) - K/dy;
            M(n,n-1) = M(n,n-1) + phi*rho*dx*dy/dt;
        end
        if j < ny
            A(n,n+1) = A(n,n+1) - K/dy;
            M(n,n+1) = M(n,n+1) + phi*rho*dx*dy/dt;
        end
        A(n,n) = A(n,n) + (K/dx+K/dy);
        M(n,n) = M(n,n) + phi*rho*dx*dy/dt;
    end
end

% 时间循环计算
for k = 1:nt
    % 计算流量
    for i = 1:nx
        for j = 1:ny
            n = (i-1)*ny+j;
            if i == 1
                q(n) = -K/dx*(h(n+ny)-h(n));
            elseif i == nx
                q(n) = -K/dx*(h(n)-h(n-ny));
            else
                q(n) = -K/dx*(h(n+ny)-h(n)) - K/dx*(h(n)-h(n-ny));
            end
            if j > 1
                q(n) = q(n) - K/dy*(h(n)-h(n-1));
            end
            if j < ny
                q(n) = q(n) - K/dy*(h(n+1)-h(n));
            end
        end
    end
    
    % 计算水头
    h = h + M\(q*dt - A*h*dt);
    
    % 可视化结果
    figure(1)
    clf
    subplot(121)
    surf(reshape(h,ny,nx)','EdgeColor','none')
    title(sprintf('Water head at t=%g',k*dt))
    xlabel('x')
    ylabel('y')
    zlabel('h')
    view(2)
    colorbar
    subplot(122)
    quiver(reshape(q,ny,nx)',0.5)
    title(sprintf('Flux at t=%g',k*dt))
    xlabel('x')
    ylabel('y')
    axis equal
    drawnow
end

上述代码中，通过定义模拟区域的大小和网格数，以及材料参数（渗透率、孔隙度、密度、重力加速度等），建立了有限元模型，并在时间循环中通过求解线性方程组来计算每一时刻的水头和流量分布。最后通过可视化结果，可以观察到地下水的动态变化过程。