matlab圆柱绕流

圆柱绕流是流体力学中的一个经典问题，可以用来研究流体在圆柱周围的流动特性。在MATLAB中，可以通过求解Navier-Stokes方程来模拟圆柱绕流。以下是一个简单的示例代码：

% 圆柱绕流模拟
clear all; close all; clc;

% 定义参数
L = 2; % 长度
H = 1; % 高度
nx = 101; % x方向网格数
ny = 51; % y方向网格数
Re = 100; % 雷诺数
U = 1; % 入口速度

% 定义网格
x = linspace(0, L, nx);
y = linspace(0, H, ny);
dx = x(2) - x(1);
dy = y(2) - y(1);
[X, Y] = meshgrid(x, y);

% 初始化速度场和压力场
u = zeros(ny, nx);
v = zeros(ny, nx);
p = zeros(ny, nx);

% 迭代求解
for n = 1:1000
    % 计算速度场和压力场
    [u, v, p] = navier_stokes(u, v, p, dx, dy, Re, U);
    
    % 绘制流线图
    if mod(n, 50) == 0
        figure;
        streamslice(X, Y, u, v);
        axis equal;
        title(['n = ', num2str(n)]);
    end
end

% Navier-Stokes方程求解函数
function [u, v, p] = navier_stokes(u, v, p, dx, dy, Re, U)
    % 计算时间步长
    dt = 0.01 * min(dx, dy)^2 / Re;
    
    % 计算中心差分
    uxx = (circshift(u, [0 -1]) - 2*u + circshift(u, [0 1])) / dx^2;
    uyy = (circshift(u, [-1 0]) - 2*u + circshift(u, [1 0])) / dy^2;
    vxx = (circshift(v, [0 -1]) - 2*v + circshift(v, [0 1])) / dx^2;
    vyy = (circshift(v, [-1 0]) - 2*v + circshift(v, [1 0])) / dy^2;
    
    % 计算对流项
    uux = u .* (circshift(u, [0 -1]) - circshift(u, [0 1])) / (2*dx);
    vuy = v .* (circshift(u, [-1 0]) - circshift(u, [1 0])) / (2*dy);
    uvx = u .* (circshift(v, [-1 0]) - circshift(v, [1 0])) / (2*dy);
    vvy = v .* (circshift(v, [0 -1]) - circshift(v, [0 1])) / (2*dx);
    
    % 计算右侧项
    Fx = -p(2:end-1,:) + Re*(uxx(2:end-1,:) + uyy(2:end-1,:)) - uux(2:end-1,:) - vuy(2:end-1,:);
    Fy = -p(:,2:end-1) + Re*(vxx(:,2:end-1) + vyy(:,2:end-1)) - uvx(:,2:end-1) - vvy(:,2:end-1) + U^2;
    
    % 计算压力场
    p_new = zeros(size(p));
    for i = 1:1000
        p_old = p_new;
        p_new(2:end-1,2:end-1) = (p_old(3:end,2:end-1) + p_old(1:end-2,2:end-1) + p_old(2:end-1,3:end) + p_old(2:end-1,1:end-2) - dx^2/dt*Fy(2:end-1,2:end-1)) / 4;
        p_new(1,:) = p_new(2,:);
        p_new(end,:) = p_new(end-1,:);
        p_new(:,1) = p_new(:,2);
        p_new(:,end) = p_new(:,end-1);
        if max(abs(p_new(:) - p_old(:))) < 1e-6
            break;
        end
    end
    
    % 更新速度场
    u(2:end-1,2:end-1) = u(2:end-1,2:end-1) + dt/Re*(uxx(2:end-1,:) + uyy(2:end-1,:) - p_new(3:end,:) + p_new(1:end-2,:)) - dt*u(2:end-1,2:end-1).*(circshift(u, [0 -1]) - circshift(u, [0 1]))/(2*dx) - dt*v(2:end-1,2:end-1).*(circshift(u, [-1 0]) - circshift(u, [1 0]))/(2*dy);
    v(2:end-1,2:end-1) = v(2:end-1,2:end-1) + dt/Re*(vxx(:,2:end-1) + vyy(:,2:end-1) - p_new(:,3:end) + p_new(:,1:end-2)) - dt*u(2:end-1,2:end-1).*(circshift(v, [0 -1]) - circshift(v, [0 1]))/(2*dx) - dt*v(2:end-1,2:end-1).*(circshift(v, [-1 0]) - circshift(v, [1 0]))/(2*dy);
end